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This paper deals with fully-connected mean-field models of quantum spins with p-body ferromag- 
netic interactions and a transverse field. For p = 2 this corresponds to the quantum Curie- Weiss 
model (a special case of the Lipkin-Meshkov-Glick model) which exhibits a second-order phase tran- 
sition, while for p > 2 the transition is first order. We provide a refined analytical description both 
of the static and of the dynamic properties of these models. In particular we obtain analytically the 
exponential rate of decay of the gap at the first-order transition. We also study the slow annealing 
from the pure transverse field to the pure ferromagnet (and vice versa) and discuss the effect of 
. . , the first-order transition and of the spinodal limit of metastability on the residual excitation energy, 

04 ■ both for finite and exponentially divergent annealing times. In the quantum computation perspective 

this quantity would assess the efficiency of the quantum adiabatic procedure as an approximation 
^3 , algorithm. 

(N 

PACS numbers: 

lO ■ I. INTRODUCTION 

Finding the minimum of a cost function defined on a discrete configuration space is the central task of combinatorial 
i-G optimization. Depending on the problem considered (i.e. the shape of the cost function), there exists or not fast 

Jn ' (running in polynomial time with respect to the number of variables) algorithms for classical computers that performs 
^ , the minimization, as classified by the computational complexity theory [l|. 

In more physical terms this problem corresponds to finding the groundstate of an Hamiltonian (cost function) 
depending on discrete degrees of freedom (spins). This analogy has suggested an optimization algorithm named 
simulated annealing [2|], that proceeds through a stochastic exploration of the phase space, according to transition 



'^ I rules obeying the detailed balance condition for a positive temperature, which is slowly reduced from a very high 
"*^ ■ value down to zero. In this way energy barriers can be jumped over through thermal fluctuations, the flnal state of 

a' the system is the equilibrium at zero temperature, hence concentrated in the sought-for minimum of the cost function. 
Provided with a quantum computer, that is a device that obeys the laws of quantum mechanics at the level of its 
1-^ ■ computing units, one can follow a similar idea but with quantum fluctuations replacing the thermal ones; this strategy 
^ , is known as quantum annealing |3i, 4j], or quantum adiabatic algorithm [5|, see jSiTi for reviews. The control parameter 
O ' that replaces the temperature allows to tune the relative strength of the potential energy (the cost function) and of 
CJ ] the "kinetic energy" (for instance a transverse field for spins 1/2). The system is initially prepared in the groundstate 
of the latter, then evolves according to Schrodinger's equation with an Hamiltonian that slowly interpolates between 
^sj ' the kinetic and the potential energy. If this interpolation is sufficiently slow the system remains at all times in the 
^ instantaneous groundstate of the Hamiltonian, and in particular at the end of the evolution it is found in the desired 
C^ minimum of the cost function. 

^^ , To assess the efficiency of these algorithms one has to specify how slow the evolution of the control parameter has to 

be in order that the final state indeed corresponds to the groundstate. In the quantum setting, which shall be the focus 
of this article, this condition is provided by the quantum adiabatic theorem [8| which, roughly speaking, states that 
f^ ' the interpolation time has to be larger that the inverse square of the minimal energy gap between the instantaneous 
groundstate and the first excited state encountered along the interpolation. It thus appears that quantum phase 
transitions Q, where the gap closes in the thermodynamic limit, constitute the bottleneck for the efficiency of the 
quantum annealing. First-order phase transitions, at which the gap is typically exponentially small in the system size, 
are in this respect worse than second-order transitions for which, at least in non-disordered systems, the gap is only 
polynomially small. 
; ^ ] Random instances of combinatorial optimization problems provide useful benchmark ensembles of cost functions 

Cd ' on which to test various algorithms [i3| ■ They have been the object of an intense research activity at the crossroad 
between computer science, mathematics and theoretical physics |11| . Several phase transitions have been unveiled 
that affect the typical number and organization of their ground and excited states p^ - ll4l | . More recently the tools 
that allowed the description of these transitions have been extended to take into account the additional effect of a 
transverse field on the corrugated random cost function [l^ \1^ . First order phase transitions as a function of the 
interpolating transverse field have been observed in some models [17, 18] ; this did not come as a surp rise as it is a 
recurrent feature of mean- field quantum disordered systems that have been extensively studied |19l - l22| . 
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Even if a first order transition in a given model means that the corresponding combinatorial optimization problem 
will only be solved exactly in a time exponentially large in the system size, many questions remain open at this 
point. First, one should try to compute the exponential rate of growth of the adiabatic time. Second, and maybe 
more importantly, one should investigate what is the final energy of an evolution that is too fast to respect the 
adiabaticity criterion (a question reminiscent of the Kibble-Zurek mechanism, see ^23.] for a recent review). Besides its 
intrisic physical relevance, this point is also deeply related to important issues in computational complexity theory, 
namely hardness of approximation results [2Jj. Indeed, for some combinatorial optimization problems (MAX-3-SAT 
for instance, or even MAX-3-X0RSAT whose decision version is in P) it is not only difficult to compute the exact 
value of the minimum cost function, but even providing an approximate answer that is asymptotically more precise 
than taking the value of the function at a random point in the configuration space is also a difficult problem j25| . 
Hence a fast non-adiabatic evolution has a computational interest if one can find a good compromise between the 
evolution time and the residual energy. In the classical case hardness of approximation results are often obtained 
via the PCP theorem |2a]. In the quantum complexity litterature a quantum analog of the PCP theorem has been 
conjectured in 12711 . For recent works on the approximation algorithms in the quantum complexity setting we refer 
the readers to [lalli]. 

In this paper we shall investigate the annealing on non-adiabatic timescales for a class of ferromagnetic, non 
disordered, mean-field models of the fully-connected type, with p-spin interactions. These can of course not be 
considered as difficult optimization problems. However, despite their simplicity that allows for an analytical resolution, 
they exhibit some of the features expected also in more realistic optimization problems, and therefore constitute useful 
toy-models to study. The statics [30l - l36j and the dynamics, both for quantum annealings |37h40| and for quantum 
quenches |4ll - l43l |. of this kind of models have been largely studied. From a technical point of view these models 
are relatively simple because their mean-field character allows for a semi-classical treatment, the small parameter in 
this limit being the inverse of the size of the system (instead of h in usual semi-classical computations) . In most of 
previous works this semi-classical limit has been achieved through the introduction of spin coherent-states |33l . |4J] , or 
instantonic computations [3J|. Here our treatment will have more of a WKB flavour, with the magnetization playing 
the role of a particle coordinate. Moreover most of these works dealed with second-order phase transitions, at the 
exception of [34. 35]. which studied the statics of models with first-order transitions. The annealing dynamics of such 
models with first-order phase transition was not investigated before, to the best of our knowledge. 

Let us now explain the structure of the paper. In Sec. |lT] we introduce the definition of the models (jll Ap and 



present their thermodynamic behavior and phase diagram (|IIBp . Section iHIl contains our investigations on their static 
properties, at a refined level with respect to the thermodynamic quantities. As explained in a first part (Sec. IIII A|) of 
this section, their mean-field character induces strong symmetries that allow to decompose their spectrum in various 
disconnected sectors. We provide an ordering theorem between different sectors in Sec. IIII BJ In Sec. IIII O we discuss 
the qualitative features of the spectrum of the model, and point to the following parts of the text where they are 
quantitatively derived. The main technical result is established in Sec. IIIIDl where we show how to determine the 
eigenvectors, at the leading level in the thermodynamic limit. This is then applied to the computation of various 
quantities: the density of states inside one sector (Sec. IIIIEp . the finite gap between levels away from transitions 
f Sec. IIII F)) . and the exponentially small gaps in Sec. IIII Gl The latter part is divided according to the location of the 
quasi-degenerate levels in the spectrum, we consider in particular the exponentially small gap between the groundstate 
and the first excited state at a first-order transition in Sec. IIII G 11 and the exponentially small gap between the two 
ferromagnetic phases for even p in Sec. IIII G 21 The dynamics of the models is studied in Sec. IIVI After a precise 
definition of the annealing procedure in Sec. IIVAI we recall the basic mechanism of the Landau-Zener model in 
Sec. IIV Bl and discuss the behavior of the dynamics that it suggests, in view of the properties of the spectrum derived 
previously. The actual results are presented in Sec. IIVDI (resp. Sec. lIVEp for annealing on exponentially large (resp. 
finite) timescales. A simplified model, introduced in Sec. IIV CI is also studied for comparison. We finally draw our 
conclusions in Sec. |Vl Some technical details are deferred to a series of Appendices. 

II. DEFINITION AND THERMODYNAMIC PROPERTIES OF THE MODELS 

A. Definition 

We shall consider Hamiltonians of interacting spins 1/2, acting on the Hilbert space spanned by 
{|2l)|o:= (cti, . . . ,crAr) e { — 1,4-1}^}. We denote af , (rf and af the Pauli matrices acting on the i-th spin, and 
recall that in this basis, fff |a;) — cri\q}, trf |ct) = la;^'-'), where ct*^*^ is the configuration obtained from a by flipping the 



i-th spin. The transverse and longitudinal magnetization per spin operators are defined as follows : 



N , N 



The Hamiltonian of the fully-connected p-spin ferroniagnet is usually defined as —N(jh^Y ~ TA^m^, i.e. with p-body 
interactions along the z axis, and a transverse field F along the x axis. The dependency in N is chosen to ensure 
the extensivity of the model in the thermodynamic limit. For future convenience we shall trade F for a parameter 
s £ [0, 1] and define 

H{s) = -Ns{rffY - N(\ - s)to^ . (2) 

Up to a change of the energy scale these two definitions are equivalent, with the correspondance F = —^^. The two 
limits s = and s = 1 corresponds to a pure transverse field and pure ferromagnetic interactions along z, respectively. 
The mean-field character of the model arises from the form of the interacting term, which depends on the total 
magnetization only. The p = 2 case corresponds to the quantum Curie- Weiss model, which can also be viewed as the 
anisotropic version of the Lipkin-Meshkov-Glick (LMG) model 30, 33, 4^ (the general LMG model contains pair-wise 
interactions in the y and z directions). The case p > 3 was investigated in [3J], and generalized in [3a] to a model 
where both m^ and fh^ are raised to arbitrary powers, and in [36| with the addition of antiferromagnetic pairwise 
interactions. The methods and results developed in this paper for the model of Eq. ^ are easily extended to these 
generalizations, as sketched in the conclusions. 

B. Thermodynamic properties 

We shall first briefiy explain how to compute the free-energy density of this model, in the thermodynamic limit, 
and discuss its phase diagram. Similar derivations can be found in jig . |3J, |36| . For a rigorous treatment of such 
models we refer the reader to [46[. The partition function at inverse temperature /3 can be obtained by mapping the 
quantum problem to a classical one with one additional imaginary time direction. Using the Suzuki- Trotter formula 
to disentangle the two non-commuting terms in the Hamiltonian, and inserting representations of the identity between 
each of the Ng Suzuki- Trotter slices one indeed obtains: 



Z(/3,s) = Tre-^^(^)= lim V ff (a(a)|e^^^(™^)''e^(i-^)^™na(a -H 1)) . (3) 



In this expression ct(1), . . . ,a{Ns) are Ng Ising spin configurations, with periodic boundary conditions a{Ns + l) = o^(l) 
As m^ is diagonal in the basis chosen one obtains 



Z(/3,s)= lim y TT e^^^(^^"- -'("))" (a(a)|e^(i-^)^"na(a + l)). (4) 

o:(l),...,<T(JV,)a=l 



Thanks to the mean-field character of the model one can reduce the problem to a single-site one by defining 
m[a) = jj'Ylii=i'^i{'^) ^nd imposing this definition, for each a, by an exponential representation of the Dirac dis- 
tribution with conjugate parameter A(a): 



f -1-4 dr7i(a)dA(a) 



AT -^E(sm(a)P-A(a)m(a))+lnTr]J - 



g„^(A(a)S- + (l-s)S-) 



(5) 



Making the natural assumption that the dominant contribution comes from values of m{Q) and A(a) that are constant 
in imaginary time and equal to to, A respectively, and evaluating the integral via the saddle-point method yields 



/(/?,s)= lim --LinZ(/3,s) = infext 

N^oo pl\ m X 



1 



-smP + Xm- -ln2 cosh(^V^2 + (i - s)2) 



(6) 
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FIG. 1: Thermodynamic properties of the p = 2 model. Left panel: phase diagram, the line indicates the value Tc(s) of the 
critical temperature for the second-order transition between the paramagnetic and the ferromagnetic phases. Center panel: 
longitudinal {rriz, red lines) and transverse {rrix, black lines) magnetizations as a function of s, for T — 1 and T = 0; in 
the ferromagnetic phase rrix is independent of T, as apparent from the last expression of Eq. ^ with p = 2. Right panel: 
free-energy density as a function of s for T = 1, and ground-state energy as a function of s; the arrows indicate the values Sc(T) 
of the transition between the paramagnet (black lines) and the ferromagnet (red lines). 



The stationarity conditions for this function of (rn, A) are 



X ^ p sm 



p-i 



A 



x/A2 + (1 - sy 



: tanh(/3v'A2 + (1 - s)2) 



(7) 



Note that an alternative derivation of this result consists in making a mean-field approximation {fh^y — )■ {fh^)P + 
p{fh^')P~^{fh^ — {m^)) in the Hamiltonian and setting self-consistently the average value in the single-spin problem 
thus obtained [43, 14^ ■ 

The various observables can be expressed in terms of the relevant critical point (?n*(/3, s), A*(/?, s)), in particular 
the longitudinal and transverse magnetization per spin read respectively 



(m^) = m^.{l3,s) , (m^) 



1 



VA.(/3,s)2 + (l-s)- 



:tanh(/3v/A*(/3,s)2 



1 -si 



-m,{P,sf-P . (8) 
P 



These expressions are easily obtained by adding to the Hamiltonian appropriate fields conjugated to the observables 
and by deriving the variational free-energy with respect to these additional fields; the last expression of Eq. ([5]) is 
only valid under the assumption m*(/3, s) ^ 0. 

Let us first discuss the solution of these equations for p = 2, and present the associated phase diagram. The point 
(to. A) = (0, 0) is always a solution of Eq. ([7]). There is however a line in the (s, (3) plane separating a paramagnetic 
phase (at low values of /3, s, i.e. high values of the temperature and transverse field) where it corresponds to the 
global minimum of the function in ([6l), from a ferromagnetic phase where it becomes a local maximum. In the latter 
phase there appears two global minima related by the symmetry operation (to. A) — >■ (—to, —A). The spontaneous 
longitudinal magnetization m^ (/3, s) > grows continuously from at the border of the ferromagnetic phase, with 
the usual mean-field exponent /3 — 1/2. The phase transition is thus of second order, the free-energy and its first 
derivatives being continuous at the transition. These properties are illustrated in Fig. [T] 

Consider now the case p > 3. The paramagnetic solution (to, A) = (0, 0) of Eq. ([7]) is then a local minimum (with 
respect to to) of the function in ^ for all values of (/3, s). For low values of /?, s this is the only minimum of ^. 
Beyond a line /3sp(s) (or equivalently Sgp (/?)), another local minimum appears discontinuously in to*(/3, s) > (if p > 4 
is even there is also a symmetric one in —m^,(/3, s)). At its appearance this non-trivial local minimum corresponds to 
an higher free-energy density than the paramagnetic one. It is only for strictly larger values of /?, s than their free- 
energy density becomes equal, on the line /3c (s) > /3sp(s) (or Sc(/3) > Ssp (/?)). The model thus exhibits a first-order 
phase transition along the line /3c (s), associated to a discontinuity in the first derivatives of the free-energy density, 
which implies in particular a discontinuity of the magnetizations. The line /3sp(s) is the spinodal of the ferromagnetic 
phase, i.e. the limit of its existence as a metastable local minimum of the free-energy. Note that the paramagnetic 
phase is always locally stable, there is thus no spinodal line for this phase. The features of this first-order transition 
are illustrated on Fig. [51 to anticipate the discussion of the rest of the paper the results displayed there are at zero 
temperature, yet they would be qualitatively identical at any positive temperature below the transition temperature 
of the classical model (at s = 1). 

In the remaining of this section we shall collect for future use some more explicit formulas valid in the zero- 
temperature limit, which will be the most useful case in the following of the paper. The groundstate energy density 
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FIG. 2: Thermodynamic properties of the p = 3 model (all cases with p > 3 are qualitatively identical). Left panel: phase 
diagram, the solid (black) line stands for the first-order transition line Tc{s), the dashed (red) line being the spinodal curve 
T'sp(s) for the limit of existence of the ferromagnetic phase. Center panel: longitudinal (m^, red line) and transverse {rrix, black 
line) magnetizations as a function of s, at T = 0. Solid thick part of the curves correspond to the thermodynamically relevant 
phase (paramagnetic for s < Sc, ferromagnetic for s > Sc), dashed light ones to the metastable ones (s > Sc for the metastable 
paramagnet, s € [ssp, Sc] for the metastable ferromagnet). Note the square-root singularity at s^p for the magnetizations of the 
ferromagnetic phase. Right panel: groundstate energy density as a function of s; solid thick and dashed light have the same 
meaning as in the center panel. 



is obtained from ^ as 



egs(s) = lim /(/3,s) 



,3-i-oo 



inf ext 

m A 



~smP + Xm- v^A2 + (1 - s)^ 



(9) 



One can solve explicitly the stationarity condition virith respect to A, vifhich yields A = (1 — s)m/\/l — m? and thus 

egs(s) 



inf 



iP~{l-s)^ 



(10) 



The energy corresponding to the paramagnetic state ?7i = is ep,„(s) = — (1 — s). For p > 3 the ferromag- 
netic phase exists when s G [ssp,l], where Ssp is the zero-temperature limit of the spinodal line. We shall denote 
771* (s) > the non-trivial solution of the stationarity equations corresponding to a local minimum for s G [ssp, 1], 
and efm(s) = —sm^,{s)'P — (1 — s)-\/l — m^,{s)'^ the corresponding energy. We also define 77ii(s) and ei(s) as the mag- 
netization and energy of the local maximum (unstable phase of intermediate magnetization) of the function in (|10l) . 
These magnetizations are the solutions, ordered with < mi{s) < m^{s), of the equation 



m = p 



1-s 



^P-i 



VT 



The average magnetizations are then given by 



(m^) = v/l-"i*(s)2 



(11) 



(12) 



At the spinodal point we call TOgp — m*(ssp) = "T-i(ssp) and egp = efm(ssp) — ei(ssp) the longitudinal magnetization 
and the energy, while the first-order transition happens for Sc such that epin(sc) = efin(sc) = ec, with rric = m.^,{sc)■ 
Explicit formulas can be given for these quantities. Consider first the spinodal point. 771^.(3) is the solution of an 
implicit equation of the form m = g{m,s), with the function g defined by the r.h.s. of Eq. (jlip . The spinodal 
corresponds to a bifurcation of this implicit equation, in consequence (rnsp,Ssp) are solutions of nigp — g(?77,sp, Ssp) 

. Solving this system yields 



and 1 = 1^ 



'sp - \/i- ^^— Y 



^sp 



1 I „ (p-2)(P-^)/^ 

-^ ~^ P(p-l)(p-l)/2 



(13) 



At the first-order transition point, the equation rric = g{'mc,Sc) is supplemented by the condition epni(sc) = efin(sc), 
which leads to 
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FIG. 3: Spinodal line Ssp(p) and phase transition line Sc(p) in the (s,p) plane, in logarithmic scale on the p axis. For p large, 
Sc — >■ 1/2 whereas Sap — > (the domain of existence of the metastable ferromagnetic phase grows with growing p). 



The dependency on p of the spinodal and critical point parameters Ssp and Sc are plotted on Fig. [31 From the above 
explicit expressions one can in particular work out the large p asymptotics, that read 






m. = l-^ + 0(p-), ^c = i-^ + Ob-) 



(15) 



III. DETAILED DESCRIPTION OF THE SPECTRUM 

We shall now turn to a refined description of the statics of the model, beyond the computation of the thermodynamic 
limit of its free-energy density. This detailed study of the eigenvalues and eigenvectors of H{s) will be crucial for the 
understanding of the annealing dynamics presented in Sec. IIVI This section is organized as follows. In Sec. IIII Al we 
exploit the symmetries of the model to decompose its Hilbert space into several disconnected sectors. In Sec. IIIIBI 
we prove a result on the relative ordering of the eigenvalues between different symmetry sectors, and provide a finer 
conjecture motivated by numerical evidences. The qualitative features of the spectrum inside one symmetry sector 
are discussed in Sec. IIII C[ and the rest of the section is devoted to the quantitative derivation of these properties. In 
Sec. IIIIDI we obtain the solution of the eigenvalue equation inside one sector, at the leading exponential level (in a 
semi-classical fashion) . This main technical result is then exploited to obtain the density of states inside each sector 
(in Sec. IIIIE| . the finite gaps between eigenvalues (in Sec. IIII ¥\ and the exponentially small gaps (in Sec. IIII Gl) . in 
particular at the first order transition of models with p > i (see Sec. IIIIG ip . and in the ferromagnetic phases for 
evenp (cf. Sec. IIII G2p . 



A. Decomposition of the Hilbert space in spin sectors 



The diagonalization, be it numerical or analytical, of a quantum Hamiltonian is in general a very difficult task 
because of the exponential growth of the dimension of the Hilbert space with the size of the system. For the fully- 
connected mean-field models under study this difficulty is greatly reduced thanks to their highly symmetric structure: 
as a matter of fact the Hamiltonian is invariant under the permutation of any pair of spin indices. 

Let us first briefly explain how to exploit this symmetry in an abstract and general way. The Hilbert space % 
of an iV-component system is the tensorial product T-L = V"^^ of the space V of each component. The theory of 
representation f49| asserts that such a tensor product can be decomposed as a direct sum of vector spaces, classified 
according to their symmetry properties with respect to permutations. More precisely, in order to construct V®^ one 
has to sum over the Young diagrams with N boxes and no more than d rows, where d is the dimension of V . Each 
of these diagrams gives rise to a Young symmetrizer, i.e. an operator on V®^ that, roughly speaking, completely 
symmetrizes along each row and antisymmetrizes along each column of the diagram. The tensor product V®^ can 
then be written as the direct sum of the images of the Young symmetrizers; the degeneracies in this sum, as well as 
the dimensions of these images, can be computed from the shape of the diagram. This decomposition can be useful 
only if the Hamiltonian itself respect such permutation symmetries, as it becomes block-diagonal once written in this 
basis. 

This general theory that we only sketched above greatly simplifies in our case, and its consequences can be under- 
stood with more physical arguments. The important point is that the dimension d of the base space of a spin 1/2 is 



only 2 here. In consequence the Young diagrams have at most two rows, and the sum over the diagrams reduces to 
a sum over the number K = 0,1, . . . , [-yj of elements in the second row. This number counts the pair of spins over 
which the antisymmetrization procedure is accomplished. It can be given a more intuitive interpretation as follows. 
The operators S" = -y™" with a = x,y, z obey the commutation rules of an angular momentum; the total spin 
operator S"^ = (S^)^ + (S*^)^ + (5*^)^ has thus eigenvalues of the form S{S + 1) with 5* integer or half-integer. It turns 
out that the images of a Young symmetrizer with a given value of K are eigenspaces of S'^, with total spin S — ^ — K . 
In particular the states of maximal spin N/2 correspond to fully symmetric states. More generally, the results of the 
abstract construction can be recovered by using recursively the standard rules for the addition of angular momenta. 
Let us now summarize these results and write explicit formulas for the matrix elements of the Hamiltonian. There 
are 

,,N fN\N+l-2K fN\ f N \ 

distinct eigenspaces of S'^ with spin S = N/2 — K (as could be expected the fully symmetric space i^ = is unique). 
Each of them has dimension 25'-|-l = iV-|-l — 2K; using the second expression of J\fJ^ an easy computation allows 
to check that the total dimension of the Hilbert space is indeed 

LfJ 

Y,J^KiN + l-2K) = 2^ . (17) 

K=0 

The Hamiltonian H{s) is stable with respect to this decomposition; moreover its action on one of these subspaces 
depends only on the value of K, not on the choice of one of the Af^ degenerate sectors. We shall denote H^-^^s) 
the restriction of H{s) to one of the subspaces of spin N/2 — K, or equivalently view H'^^\s) as a square matrix of 
order iV + 1 — 2K. We will also use the notation k — K/N . The subspace on which H^-^\s) acts is spanned by the 
basis of eigenvectors of rh^, written \m; K)^ with the iV + 1 — 2K possible values of m,: M^ = {— 1 + 2k, — 1 + 2fc + 
2/N, . . . , 1 — 2fc — 2/N, 1 — 2fc}. The action of ffv^ on this vector amounts to increase or decrease the value of m by 
its minimal amount 2/N , i.e. 

^{m;K\Trv'\m';K)^^ = - ^(1 - 2fc + max(m,m'))(l - 2fc - min(m, m')) for |m - m'| = 2/iV . (18) 

The matrix representing H^^'{s) in this basis has thus a symmetric tridiagonal form, with matrix elements: 

^{m;K\H^'^'^{s)\m;K), = -N smP , (19) 

^{m;K\H^^\s)\m';K)^ = -ivi— ^ ^(1 - 2A: + max(TO, m'))(l -2k- min(m, m')) for |to - m'| = 2/A^ .(20) 

One can also define a second basis spanned by the eigenvectors \m; K)x of ffv^ . The expression of m^ in this basis is 
nothing but Eq. (jlSp with the interversion of the indices z and x. In this basis the matrix representation of H^^' (s) has 
a diagonal part corresponding to the action of the transverse field; the interaction term {fh^Y has a band diagonal form, 
with non-zero matrix elements between eigenvectors \m; K)x and |to'; K)x when ^-{m—m') G {p,p—2, . . . , —p+2, —p}. 
We shall give and use their explicit form, in the large N limit, in Sec. IIIIDI 

These symmetry considerations thus allow to reduce the complexity of the full diagonalization of the Hamiltonian 
from a matrix problem of size 2^ to [-jj + 1 matrices of sizes at most N + 1. This great simplification will be used 
in the following both for numerical and analytical computations. 

The above reduction is valid for any model symmetric under all permutations of spins, irrespectively of the precise 
form of the interactions. The Hamiltonian of Eq. ([2]) exhibit additional symmetries: 

• for odd values of p the spectrum is invariant under the transformation E — > —E. Consider indeed 
an eigenvector |-0) of iJ'-^'(s), with eigenvalue E, written as \ip) — ^ c„i\'m;K)z- Then the vector 

\ip') = ^ c-„ii—l)~"^\m; K) z is an eigenvector of H'^^\s), with eigenvalue —E. 

• for even values of p the Hamiltonian is symmetric under global longitudinal magnetization reversal, and this 
implies that each iJ ^^^ (s) can be further decomposed in a block diagonal form, with two blocks of sizes [ ^J + 
1 — K and [^1 ~ ^- O^^^ can justify this statement in two ways. When acting on a basis vector \'m;K)x, 
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FIG. 4: The eigenvalues -E{^'(s) of the Hamiltonian H{s), for TV = 12, p = 2 (left) and A^ = 12, p = 3 (right), 
values of K are distinguished by different colors and line styles. 
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the operator (rh^)P, with p even, produces a vector whose decomposition over \m';K)x is non-zero only for 
magnetizations m' such that N/2{m — m') is even: non-zero off-diagonal matrix elements in the x basis are only 
found at an even distance from the main diagonal, hence in that basis the parity of the number of spin flips 
with respect to the fully polarized vector in the x direction is conserved by the Hamiltonian. 

In addition, the matrix representing H^^'{s) in the \m]K)z basis commutes with the matrix with 1 on the 
anti-diagonal (from bottom left to top right), that represents the reversal of the magnetization along the z axis. 
The eigenvectors of H'^^^ (s) can thus be divided between those that are symmetric or antisymmetric under this 
transformation. 



B. Ordering properties of the spectra 



For each value of K the restriction H'^^^ (s) of H{s) to a subspace of spin N/2 — K has A^ + 1 — 2K real eigenvalues. 

We leave implicit the dependency of these quantities on 



{K) 



?{K) 



< Efl,A^). 



that we shall denote E^'^^s) < EY^'{s) < 

N and p, that are understood to be fixed in this whole subsection. By construction the Hamiltonian has no matrix 

elements between sectors of the Hilbert space corresponding to different values of K; one could thus a priori think 



AK) 



?iK 



that the spectrum {E^ (s)} has no relationship to {E^ (s)} for K ^ K' . A quick look at the numerical results 
displayed in Fig. 2] reveals on the contrary that there are strong ordering rules between the energy levels of different 
spin sectors, reminiscent of the Lieb-Mattis theorem for the antiferromagnetic Heisenberg model [SOl (see also [51[ for 
a more recent treatment of the ferromagnetic case) . As a first step we shall prove that the groundstates of each sector 
are strictly ordered according to the spin of the sector, i.e. that 



.(0)^ 



.(1) 



(LfJ)/ 



E^,'>{s)<E^^\s)<...E'^^^>{s) V.e[0,l] 



(21) 



in such a way that the global groundstate of ii{s) lies in the fully symmetric subspace K = Q, oi maximal spin N/2. 
The proof goes as follows. Let us denote Jt/i) the eigenvector of H^^^ corresponding to its groundstate eigenvalue 
Eq , for some value oi K > Q. We decompose this vector on the basis in which m^ is diagonal. 



y^ Cm\m;K)^ 



(22) 



meM'S. 



We saw above that in this basis i/^^' is a tri-diagonal matrix whose off-diagonal elements are all positive (cf. Eq. ([20)) ). 
The Perron-Frobenius theorem thus ensures that the coefficients Cm can be chosen to be all strictly positive. Let us 
now define a vector |^') belonging to the space on which iy(^~i) acts, according to 



W)= Y. c™|m;ir-l) 



(23) 



meM'S. 



In a column representation this amounts to supplement {ip) with two null rows corresponding to the two values 
m = ±(1 — 2k + 2/N). This vector being normalized, the variational principle asserts that E'q < {tp'\H'^^~^^\il)'). 

One finds easily that 

meA^g\l-2fc 
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(24) 



The difference of the square roots being strictly positive, as well as the product CmCjj^jf.2., one thus obtains 

^0 i^) "^ -^0 (■*) ^s lo'^S ^^ ^ < 1- ^^ the other hand for s = 1 the matrices are diagonal and the ground- 
state is obviously Eq '{s = 1) = —N{1 — 2kY, which also obeys the strict inequality Eq ~ '{s = 1) < E^ -^(5 = 1). 
This completes the proof of Eq. ([?T|) . 

A closer look at the plots in Fig. H] suggests that not only the groundstates are ordered between one sector and 
another, but also that excited states are interleaved in a regular way. For instance the first excited state of one sector. 
El , seems to always have a lower energy than the groundstate of the following sector, Eq ' . More generally, we 
propose the following conjecture based on this numerical investigation: for all n < N/2 if p is odd, n < A^ if p is even, 
one has 

4°)(s) < Ei'l.is) <■■■< 4"'(.) V. G (0, 1) . (25) 

In particular, if this statement is true, the first excited state of the Hamiltonian H{s) in the full Hilbert space is 
always in the sector of maximal spin. The consistency of this conjecture for s close to and 1 can easily be checked 
by perturbative expansions. 

C. The salient features of the spectrum 

In this section we describe qualitatively the main features of the spectrum of eigenstates in the symmetric sector 
of maximal spin (all sectors behaving in a similar way), that are apparent by visual inspection of the figures. We 
emphasize the connections with the thermodynamic computations of Sec. IIIBi and also point to the following parts 
of the article where these properties are derived quantitatively. 

Let us begin with the p = 2 case, for which the zero-temperature limit of the thermodynamic computation predicts 
a second-order phase transition at Sc = 1/3. The complete spectrum of the symmetric sector is plotted on the left 
panel of Fig. [S] for N = 60. One observes indeed a good agreement with the shape of the groundstate energy predicted 
previously and plotted in the right panel of Fig. [1] Looking more carefully at the two states of lowest energy, one sees 
that the gap between them is of order 1 (in extensive energy E = Ne) in the paramagnetic phase (i.e. for s < Sc), 
but exponentially small in N in the ferromagnetic phase and indistinguishable in this figure. This exponentially 
small splitting is the consequence of the existence of the two magnetizations ±m^{s) minimizing the thermodynamic 
groundstate energy pH)) . On the right panel of Fig. [5] we display the gap between the lowest states for two finite values 
of N, along with the analytical prediction of Eq. ([5^ for its limit when N —> cx), that we shall obtain in Sec. IIIIEI 
The rate of the exponential splitting in the ferromagnetic phase will be derived in Sec. IIII G 2[ see Eq. (|M)) and right 
panel of Fig. 1111 The square-root vanishing of the gap when s ^ s~ and the behavior of the rate of exponential 
splitting when s — > s+ leads, with a finite-size scaling assumption explained in Sec. IIIIG 21 to a polynomial closing 
of the gap as A^"^/"^ in the critical regime s « Sc- This behaviour has been first predicted on the basis of the scaling 
analysis in [30] ; the lifting of the degeneracy between ferromagnetic states was also studied in [S^] ■ Note finally that 
the quasi-degeneracy of ferromagnetic states also occurs for excited eigenvalues: on the left panel of Fig. [5] one gets 
the impression that there are twice as less states on the right of a diagonal line e = — (1 — s) that on the left. This 
visual impression shall be confirmed in Sec. IIII G 31 

Let us now turn to the p = 3 case, which has a first-order transition at Sc- The spectrum of its maximal spin sector 
is displayed for A^ = 60 on the left panel of Fig. |6l The slope of the groundstate energy is discontinuous at Sc, as 
in the thermodynamic computation (compare with the right panel of Fig. ^ . At variance with the p = 2 case the 
gap between the groundstate and the first excited state remains of order 1 until one gets very close to Sc', this is best 
seen on the right panel of Fig. [Bl which displays a blow up of the lowest energy states around Sc- The minimal gap 
(reached in Sc{N) which goes to Sc in the large A^ limit) is indeed exponentially small at the first order transition; its 
exponential rate of decrease, which has been determined numerically and via an instantonic computation in |34| . will 
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FIG. 5: Left panel: the spectrum of the symmetric sector of the p — 2 model for A'^ — 60, obtained by numerical diagonalization. 
Right panel: the gap between the groundstate and the first excited state for p = 2. The solid line is the result of an analytical 
computation, see Eq. ()52p . the symbols have been obtained by numerical diagonalization. 



be computed in Sec. lIII G H and given as an explicit analytic formula in Eq. (|62p . The level repulsion between the two 
lowest cigenstates is at work only in a small neighborhood of their point of avoided crossing, and it is tempting to infer 
from this plot that the first excited eigenvector for s > Sc{N) is the continuation of the groundstate eigenvector of 
s < Sc{N) (and leads thus to the same thermodynamic observables), and vice versa. This pattern of avoided crossing 
looks actually very familiar, and can be observed in many examples involving the eigenvalues of an operator depending 
on an external parameter, s here. Let us recall the fundamental reason behind the universality of such a pattern, and 
the justification of the continuation intuition (a detailed discussion of this point can be found in [52|). The matrix 

elements of the operator H{s) are analytic functions of s. As a consequence the eigenvalues E^ (s) (in any symmetry 
sector K) are the roots of an algebraic (characteristic) equation of order iV + 1 — 2K, whose coefficients are analytic 

functions of s. Then one can prove (see for instance theorem XIL2 in [53|) that the E^ (s) are analytic functions of s, 
with at most algebraic branchpoint singularities at the (a priori complex) values of s where the roots are degenerate. 
An avoided crossing is thus due to the two eigenvalues being strictly equal when an infinitesimally small imaginary 
part is added to s. By performing the interpolation between s < Sc{N) and s > Sc{N) via a detour in the complex 
plane avoiding the branchpoint singularity one can thus define in a precise way the first excited eigenvector on one 
side of the avoided crossing as the analytic continuation of the groundstate on the other side. The thermodynamic 
calculations of Sec. Ill Bl suggest that the paramagnetic state of energy e = — (1 — s) is metastable for all values 
s G [sc, 1], and indeed one observes on the plots of Fig. [5] and [7] a continuation of this eigenstate across a series of 
avoided crossings. On the contrary the ferromagnetic state which corresponds to the groundstate for s > Sc only 
exists down to a spinodal point Sgp, beyond which the analytic continuation cannot be performed anymore. This is 
illustrated on the two panels of Fig. [71 The exponential rate of closing of the gaps encountered along the continuation 
of the paramagnetic and ferromagnetic state shall be computed analytically in Sec. lIIIG4l along with a determination 
of the area in the (s, e) plane where avoided crossings do occur (see Sec. IIIIG3p . 

All odd values of p > 3 yield behaviors similar to the p = 3 case. The models with an even value of p > 4 
exhibit both the first-order phenomenology of the p — 3 case and the exponentially small splitting between the two 
ferromagnetic states allowed by the spin-fiip symmetry. This shall be further discussed in Sec. lIIIG^ and Sec. IIII G3\ 



D. The semi-classical solution of the eigenvalue equation 



We shall now determine the structure of the eigenvectors of H{s), in the thermodynamic limit N — >■ (X), with a 
calculation formally similar to the WKB semi-classical treatment of quantum mechanics, the size of the system TV 
playing the role oi H~^. Similar semi-classical analysis have been performed for mean- field spin models in |33l l44l l54l| 
using a spin-coherent-state representation; at variance with these works we shall use here the eigenbasis of m^ or m^ . 

As explained above this computation amounts to diagonalize the matrices of order iV -I- 1 — 2K representing the 
restriction H^-^^s) to a sector of spin N/2 — K. For simplicity, and because this will be the most useful case in the 
following, we shall concentrate here on the fully symmetric K — Q situation. The generalization to higher values of K 
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FIG. 6: Left panel: the spectrum of the symmetric sector of the p — 3 model for A'^ — 60, obtained by numerical diagonalization. 
Right panel: a zoom on the lowest energy part of the spectrum around the avoided crossing between the lowest energy eigenstates 
at Sc{N). In the thermodynamic limit the first-order transition occurs at s = Sc. 
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FIG. 7: Left panel: a part of the spectrum for p — 3, with N = 320, in the neighborhood of the spinodal point of the 
ferromagnetic phase. Right panel: a zoom on the area of the black box of the left panel. All the crossings are actually avoided, 
but some of the gaps are too small to be distinguished on the picture. 



is straightforward and sketched in Sec. IIIIEl Let us look for an eigenvector |0(s,e)) of H^^'>{s) with eigenvalue Ne, 
written in the z-diagonal basis as 



l<^(s>e))= ^ (l){m,s,e)\m;0). 



(26) 



Using the expression of the matrix elements given in Eqs. (J19l20p . one obtains the equation obeyed by the coefficients 
(j){m, s,e): 



e(f>{m,s,e) — —sm^ (j){m, s,e) — \ 1 — m"^ -\ {1 — m) <j) I m -\ ,s,e 



—J'^ -m2 + — (1+m) 0f m- — ,s,e 



(27) 



To deal with the N -^ oo limit we shall make the following Ansatz on the behaviour of the eigenvector components: 
(j){m,s,e) = e~ ^^™'''^''^\ with ip{m,s,e) an a priori smooth complex function. Then 0(m + 2/N,s,e) + (j){rn — 
2/N,s,e) = 20(m, s,e) cosh(2(^'(m, s,e) + 0{1/N)), where ' denotes the derivative with respect to m, and (P7)) can 



12 



be rewritten, at the leading order, as: 



e = — sm^ — (1 — s)v 1 — m^ cosh(2iy9'(r7i, s, e)) . (28) 

Inverting this relation yields 

,, .1 , / e + smP \ 

cz) TO, s,e) = - argcosh , , (29 

^ ^ 2 ^ V (l-.s)Vl-TO2y ^ ' 

the analog of the eikonal equation in the semi-classical one-dimensional quantum mechanics context. Several points 
have to be precised for this equation to unambiguously determine the eigenvector rate function ip{m^s^e). First of 
all, not all values of e should correspond to an authorized eigenvalue of the Hamiltonian. Then the value of ip has to 
be fixed at one point to to reconstruct (p from its derivative. Finally argcosh is a multi- valued function, hence one 
should precise which of its branches to use. 

These ambiguities are actually solved by imposing that (j) is normalizable in the large N limit, and that (/?' is 
continuous in to (the coefficients of the eigenvalue equation ()27p being smooth in to.). By fixing for instance the norm 
of (j) to be of order 1, the first requirement imposes that mi[^ip{m,s,e)] = (we denote ^z and Qz the real and 

m 

imaginary part of z). To precise the meaning of the argcosh function let us first define the functions ach(t) as the 
reciprocal of cosh that maps the interval t G [1, oo) to [0, oo), and acos (t) the reciprocal of cos that maps t £ [—1, 1] 
to [0, tt]. Then, as the argument of the argcosh in Eq. (|29p is always real, it is enough to define 

'j7r±ach(-i) ifi<-l 
arg cosh t = •^ i acos (i) ifie[— 1,1] . (30) 

±ach (i) if t > 1 

There are two branchpoints in ±1, where the function thus defined is continuous independently of the choice of the 
sign of its real part. 

Let us now derive from these considerations the authorized value of the eigenvalue (per spin) e. Notice first that 
the argument of argcosh in Eq. (09)) diverges in to — >■ ±1, hence it cannot be confined to [—1,1] for all values of 
m G [—1,1]. There remains two cases to consider: either the argument crosses at least once one of the two branch- 
points ±1, or it remains larger (in absolute value) than 1 for all m. In the latter case one cannot change branch and 
the sign of the real part of ip' is constant on ?ti S [—1,1] (otherwise ip' is not continuous), hence one cannot fulfill the 
condition mi[^ip{m, s, e)] = in a non-trivial way. A moment of thought reveals that on the contrary in the former 

rn 

case one can construct a normalizable eigenvector. This implies that the range of authorized values for e is 



Image [—smP — (1 — s)v 1 — to^] U Image [—smP + (1 — s)v 1 — rri^] ■ (31) 

me[-l,l] me[-l,l] 

In particular the groundstate energy, for a given value of s, is obtained from this reasoning as 



egs{s)= inf -smP~{l~s)y/l-m^ , (32) 

me[— 1,1] L J 

in perfect agreement with the thermodynamic computation of Sec. IIIBi see Eq. (fTO|) . 

In the following sections llIIElHII G ll and lIIIG2l we shall show explicitly, in various cases, how to choose the correct 
branches of the argcosh function when crossing a branchpoint and how to determine the rate function (p{m, s, e) by 
integration of Eq. (j29p . A particularly important issue will be the occurence of multiple valid eigenvectors correspond- 
ing, at the leading order, to the same eigenvalue e. Before that we shall present a very simple example to check the 
above computation, and an alternative formulation in another basis. 

A simple consistency check of Eq. (P^ can be performed for s = 0, i.e. in a pure transverse field. In that case it is 
easy to see that for all N the groundstate has energy —N, with the eigenvector 

|</)(s = 0,e = -l)) = -^ ^ J(./i;™)l™;0). , (33) 

corresponding to all spins aligned in the x direction. These values of 4>(m, s — O^e — —1) solve exactly Eq. (j27l) : with 
the help of the Stirling formula one obtains the value of ip in the iV — > oo limit, 

1 + TO 1 — TO, 

ipo(m) = (p[m, s — 0,e — —I) = — - — ln(l -I- to) H — ln(l — to) . (34) 
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Let us now check that the computation presented above gives back this result. We have from Eq. ^^ 



(^'(m,s = 0,e = -1) = -arg cosh f^====j . (35) 

The argument of the argcosh reaches the branchpoint 1 only in m = 0. To enforce the condition mi[iRip{m,s,e)] — 
one has to choose the branches as 

(p'(TO,s = 0,e = -1) = sign(TO)-ach f A , (36) 

hence upon integration with the boundary condition ip{m = 0, s = 0, e = —1) = 0, 

1 r^^ / 1 \ 1 r^^ 1 /""^ / ^ -L m'\ 
iplm, s = 0, e = —1) ~ — sign (m) ach — ) dm' — — argtanh(77i')d?n' — — In ) dm' , 

2 Jo VVl-m'V 2 Jo 4 Jo \l-m'J 

(37) 
in agreement with the direct computation yielding (j34p . 

We shall finally present a similar computation of the eigenvectors of H^^'{s), but using now the x basis, namely we 
write 

\<l^{s,e)) = X! <l}x{in,s,e)\m;0)^ . (38) 

The coefficients (px obey the following equation (the equivalent of Eq. (I?7l) in the z basis): 

e(j)x{m,s,e) ^ -{I- s)m(j)xim,s,e)- si — - — j ^f j c/y^ im+ —{p - 2u),s,e\ , (39) 

in which we have dropped some irrelevant terms of order 1/A^. As above we look for a solution of this equation under 
the form ipxirn, s, e) = g-^VxlmjS.e)^ ^^^ ^^^ fhai the leading behaviour of (px is ruled by the equation 

e = -{l-s)m-s(^—^\ ^CP\e^v'Am,s.e)Y-u(^^-2v'Am,s.e)Y (^q) 

= -(l-.s)m-s(l-m2)P/2cosh(2^;(m,s,e))P . (41) 

This yields finally the equivalent of Eq. (j29p : 

^ 1 ^.(f e + {l-s)mV'^\ 

^^(m,s,e) = -argcosh I ( — ^ ^_^ ^ j I . (42) 

This equation suffers from the same kind of ambiguities as Eq. (1^^ . the p-th root and the argcosh function being 
multi-valued. However these ambiguities can also be solved with exactly the same reasoning as the one following 
Eq. (|29p . In most of the paper we shall use the z-basis computation; the use of the x-basis will however reveal useful 
in Sec. IIVE2I 

E. The computation of the density of states inside one symmetry sector 

We shall now present the first application of the above computation of the eigenvectors, that will give an explicit 
formula for the integrated density of eigenvalues (similar results for the LMG model can be found in [3j]). Let us 
first define this notion precisely, and emphasize its difference with another, maybe more usual, related concept. The 
full Hilbert space of the model ^ is 2^ dimensional, and its "density of states" can be defined as the microcanonical 
entropy a{s,e), such that e^'^^'^'^^de gives, at the leading order, the number of eigenvalues of ^ in the interval 
[iVe, N{e + de)]. This quantity is obtained from the free-energy density (jG]) via a Legendre transform between e and 
/3. In this section we shall however investigate a finer quantity, namely the density of eigenvalues for the restriction 
H'^^^ of the Hamiltonian to one symmetry sector (which has A^ -I- 1 — 2K eigenvalues). Consider for instance the 
fully-symmetric sector, and define 

2)o(s,e) = hm -^\{j\Ef\s) < Ne}\ (43) 
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as the integrated density of states inside that sector. We shall see at the end of this section that the knowledge of the 
integrated density of states of all sectors K = 0, . . . , [^J provides a much more detailed information on the system 
than the microcanonical entropy. 

There is actually a simple relation between I?o(e, s) and the leading order computation of the eigenvectors of the 
previous Subsection, based on the following observation: H^^\ expressed in the z-basis, is a symmetric tridiagonal 
matrix with all elements next to the diagonal of the same sign (negative) . This implies that the ordering in energies of 
its eigenvalues corresponds to the number of nodes of the associated eigenstates, exactly for the same reasons as the 
n-th excited eigenstate of a one-dimensional quantum particle described by the Schrodinger equation has precisely n 
zeroes. A proof for the discrete case can be adapted from the usual reasonings in the Schrodinger case (see [53 for a 
similar derivation in another context), and shows that the groundstate of H^^' is a Perron- Frobenius vector whose 
elements can be taken all positive, while its first excited state presents exactly one "domain wall" between two sets of 
values of m where the eigenvector is positive/negative, and so on and so forth. As we defined (l){ra) — e^^'^(™', the 
fact that excited eigenstates exhibit alternating signs translates into Lp acquiring an imaginary part. More precisely, 
each "domain wall" between opposite signs for (j){m) corresponds to an increase of its phase by ±7r. One can thus 
count the number of nodal points of 4> by integrating the imaginary part of (p' on to S [—1,1], and deduce from it the 
number of eigenvalues that have lower energies. This reasoning yields the following formula, 

2?o(sie) = — / dm3(y9'(m, s, e) , (44) 

T^ J-i 

as we have chosen in Eq. pop a branch of arg cosh with positive imaginary part. Using the value (|29p for the derivative 
of ip one obtains a completely explicit formula for the integrated density of states in the maximal spin sector, 

r-l 



X'o(s,e) 



1 

2^ 



dm 



acos 



(1 - s)V\ 



{l-s)^/\ 



-ttI - 



e[-i,i] 



< -1 



(45) 



where we defined 1(A) to be 1 if A is true, otherwise. Deriving this expression with respect to e one can equivalently 
obtain an expression for the density of states, 

/9o(s,e) = :^ / dm — I 



1 

2^ 



ehi,i] 



(46) 



y(i^^7)2(r^r7;T2yT^7qr7;;^ \^ (i-s)VT^^ 

We give in Fig. |S]some examples of the construction of ^p, for p — "i and s — 0.3, i.e. in the paramagnetic phase. 
For the ground-state energy e = — (1 — s) the argument of the arg cosh function in Eq. ([^^ is always > 1, with a 
single point of equality in m = 0. As a consequence the corresponding solution for if is everywhere real, see left panel 
of the figure. On the contrary for a slightly higher value of the energy (e = Cgs -I- 0.1 on the figure) the branchpoint 1 
is crossed at two values of to, hence the imaginary part of (p grows on this interval, on which the real part vanishes 
identically. We also present on Fig. IH] the curves for the integrated density of states Vq for p = 3 and two values of s, 
0.3 and 0.6 (in the latter case one observes a singularity at the crossing of the energy of the metastable paramagnetic 
phase) . The agreement with the density of states obtained by numerical diagonalization is very good already for small 
values of N {N — 40 on the figure). 

Let us briefiy mention here one application of this computation that shall be useful in the analysis of the annealing 
dynamics. From the integrated density of states Vq^SjC) one can define "iso- integrated density lines" eiso(s) by 
imposing that 2?o(s, eiso(s)) remains constant when s is varied. These lines correspond to the thermodynamic limit 
of the energy density of some (excited) eigenvalues, as long as no level crossings occurs. 

We shall now explain how to generalize the computation of the integrated density of states to sectors of arbitrary 
spin N/2 — K. The matrix of size iV -I- 1 — 2K representing the restriction H^ "> of the Hamiltonian to this sector 
has matrix elements given in Eqs. (J19l20p . Denoting k — K/N, one can look for eigenstates of H^^' under the form 
^-Nip(m.,s,e) ^ where the longitudinal magnetization to, is now restricted to [—1 + 2k, 1 — 2k], and the function tp is 
solution of a generalization of Eq. (1^^ , namely 



1 / e + smP 

LP (to, s,e) — - arg cosh = 

2 ^ V (l-s)v/(l-2fc)2 



(47) 



The rest of the computation follows strictly the reasoning made for the symmetric {K — 0) sector. In particular the 
groundstate energy density for a sector with k = K/N reads in the thermodynamic limit 



°ss \^) ^ 



inf 

me[-l+2fe,l-2fc] 



-s-mF - (1 - s)V(l-2fc)2-TO2 



(48) 
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FIG. 8: Plots of the eigenstate function ip, for p = 3 and s — 0.3. Left panel: for the groundstate energy, e = egs, (/? is purely 
real. Right panel, for a slightly larger energy, e — Bgs + 0.1, ip acquires an imaginary part. 
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FIG. 9: Plots of the integrated density of states Do for p — 3, s = 0.3 and s = 0.6. The solid lines are our analytical predictions 
from Eq. (jlSp : the symbols are the results of numerical diagonalization of systems with A'' = 40. 



and the density of states in that sector is 



ZTT 



1-2*; 



dm 



-l+2fc 



e + .s m^ 



e + s m 



p 



(l-s)^(l-2fc)2-m2 



(1 - s)^(l-2fc)2-TO2 

e + s m'P 



{i-s)^ii-2ky 



ehl,i] 



< -1 



(49) 



We shall finaly compare the amount of information contained in the densities of states I?fc(s, e) on one hand, and 
the microcanonical entropy a{s, e) on the other. The latter being the Legendre transform of the free-energy, we shall 
equivalently discuss this quantity. Using the decomposition of the Hilbert space into symmetry sectors the partition 
function can be written as 



LfJ N-2K 

K=0 j=0 



(50) 



where J\f^ gives the number of degenerate representations of spin N/2 — K (see Eq. (|16l) ). and E- ' is the j-th 

eigenvalue of H^^'. In the thermodynamic limit the degeneracy M^ grows exponentially with N\ on the contrary the 
sum over the eigenstates of one sector contains only a linear number of terms, and is thus dominated at the leading 
exponential order by the greatest of these terms, i.e. the groundstate energy of the corresponding sector. This leads 



to the following expression for the free-energy density, 



/(^, s) = inf 



4? (■5) - ^(-A:lnfc - (1 - k) ln(l - k)) 
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(51) 



A short computation based on the expression of egg given in (1481) reveals the agreement between this expression and 
the one obtained in Sec. Ill Bi f see Eq. ([S])). What we want to stress here is that the only "microscopic" (i.e. at the level 
of eigenstates) input of the computation is the energy density of the groundstate in each sector. The microcanonical 
entropy is thus entirely dominated by the effect of the degeneracy N^ of the various spin sectors, and is completely 
insensitive to their internal structure beyond their groundstate energy density. 

F. The computation of finite gaps 

One can estimate the energy gaps between eigenvalues from the density of states obtained in Eq. (pS)) : in the interval 
[e, e + de] of (intensive) energy one finds Np{s, e)de eigenvalues. Assuming these levels to be equispaced, the gap 
between two successive eigenvalues is, in extensive energy, l/p(s, e) [33|. This computation can be performed in any 
part of the energy spectrum; for simplicity we shall only state some results, obtained by combining this observation 
with the explicit expressions of the density of states P5|) and of its integrated form (H5|) , in the most relevant regions 
of the spectrum. 

For p = 2, i.e. in the Curie- Weiss model, one obtains in the paramagnetic phase (s < Sc) for the gap between the 
groundstate and the first excited state: 

lim [e['^ is) - 4°) is)] = , \ .. = 2^/3v/(l-s)(5c-s) , (52) 

as found in [30| , and plotted on the right panel of Fig. [Sj Note the square-root closing of the gap at the second-order 
transition. The same computation performed in the ferromagnetic phase (s > Sc) yields 



p(s,egs(s)) 



V3v/(l + 5)(s-Sc) . (53) 



This should however not be interpreted as the gap between the first two eigenvalues, but rather as half the gap between 
the groundstate and the second excited state. Indeed the level splitting between the two lowest states is exponentially 
small in N (as will be computed in Sec. IIIIG2P and the density of states does not distinguish them. In other words 
the hypothesis of equi-spacing of eigenvalues is strongly broken in this situation. 
Consider now the case p > 2. In the paramagnetic phase (s < Sc) one finds 

hm [£;("' (.) - 4") (.)] = . \.. - 2(1 - ,s) (54) 

JV^oo p(s,egs(s)) 

for the gap between the two lowest levels, which is the same result as would have been obtained if the Hamiltonian 
contained only the transverse field term. Note also that the gap thus computed remains positive at the first-order 
transition; the exponentially small gap (to be determined in Sec lIIIGT]) cannot be detected by the density of states. 
The computation of the density of states po at the groundstate energy can similarly be performed in the ferromagnetic 
phase (i.e. for s > Sc). For odd values of p one finds 

lim [e[°^ is) - 4"^ (s)] = , \ ,, = 2psm, (s)^-^ ^ip - l)m,is)^ - ip - 2) , (55) 

jv-^oo p(s,egs(s)) 

which is positive at Sc- For even values of p this computation, as explained above in the case p = 2, gives an 
information only on the gap between the groundstate and the second excited state. After a short computation one 
obtains a similar formula, 

hm [E^^^ is) - eI°^ is)] = . % .. = 2psm, (s)^-^ v/(p - l)m,(s)2 - (p - 2) ; (56) 

AT-i-oo p(s, egs(s)) 

one could expect to find an additional factor 2 in this expression with respect to the odd p case, however this factor 
compensates because of the contributions of the two minima in ±m^is) in the density of states Vq. 
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For p > 2 a richer behaviour is displayed in the neighborhood of the spinodal point of coordinates (ssp,esp). In 
particular in the limit s — > s~ one finds after some computations a scaling behaviour for the integrated density of 
states, of the form 

2?o(ssp - Ss, Csp + Se) - I'o(ssp, esp) ~ Ss g{{Se + ef^{ssp)6s)6s^^/^) , (57) 



where Sfyj^issp) — —"if + ^/l — m"^ is the derivative of the energy of the ferromagnetic metastable state at the 

spinodal. The scaling function g{z) is monotonously increasing, behaves as \z\^'^ for \z\ —>■ oo, and vanishes in 
one point we shall denote zq. The iso- integrated density line that goes through the spinodal point behaves thus as 
Giso{ssp~Ss) ^ ~ef^{ssp)6s + zqSs^'^ . Moreover in the scaling regime the density of states can be obtained by deriving 
the above relation, namely 

poissp - Ss, e,p + Se) ^ Ss'^/' g'{{Se + e',^^{s,p)Ss)Ss-^/^) . (58) 

In consequence the finite gap between the eigenstate level that reaches the spinodal point and the first excited state 
above it closes when s — 7> s^ as Ss^/^ /Q'{zo). This fifth root is to be contrasted with the square root singularity 
for the groundstate of p = 2. Moreover we should warn the reader that the correction zqSs^^^ in the expansion of 
the eigenstate energy is crucial: in z = the scaling function Q is finite but has no derivative, and the expansion 
of po without taking into account the correction leads to po('Ssp — Ss, egp — efjjj(ssp)(5s) ex Ss~^''^ , which modifies the 
exponent from 1/5 to 1/4. 

G. The computation of exponentially small gaps 

As discussed qualitatively in Sec. lIIICJ the energy gaps between two successive levels are, in some regions of the plane 
(s, e) depending on p, exponentially small in the size A^ of the system. This section is devoted to the computation 
of the exponential rate of closing of those gaps. It is divided in four parts; we shall first investigate the avoided 
crossing between the groundstate and the first excited state at the first-order transition of the models with p > 3 (in 
Sec lIIIGT]) . then compute the exponentially small splitting between the two lowest levels in the ferromagnetic phase of 
even p models (in Sec. lIIIG^ . The next two subsections will be devoted to exponentially small gaps between excited 
states; in Sec lIIKTBl we shall determine the values of (s, e) where these avoided crossings do occur, and in Sec. lIIIG4l 
we will concentrate on the avoided crossings encountered by the metastable continuations of the paramagnetic and 
ferromagnetic groundstates. 

From a technical point of view the common pattern behind the appearance of an exponentially small gap is the 
existence of two valid solutions of the semi-classical eigenvalue equation ((29ll for the same value of e. This approx- 
imate degeneracy is lifted at the exponential order, the splitting between the two levels being proportional to the 
exponentially small scalar product between the two quasi-eigenvectors computed at leading order. This can be shown 
by writing the eigenvalue equation in the two dimensional Hilbert space spanned by the two quasi-eigenvectors. In 
more physical terms this corresponds to the semi-classical approximation of quantum mechanics for a double-well 
potential, in which the two lowest energy levels have a gap exponentially small in \/h. 

1. The exponentially small gap at the first order transition for p > 3 

The first case we shall consider is the exponentially small gap between the groundstate and the first excited state 
at the first-order transition for p > 3. We first concentrate on the odd p case for simplicity, the modifications to be 
made when p is even are discussed afterwards. 

From the thermodynamic considerations of Sec. lIIBl we showed that this transition happens at a (p-dependent) value 
of s denoted Sc, where the infimum in the definition (llOp of the groundstate energy is reached for two distinct values 
of TO, i.e. in TO = and to = TOc > (see Eq. p4[l for the values of Sc and to.c). The paramagnetic and ferromagnetic 
phases have thus the same energy Cc. In terms of the eigenstate computation, this observation translates into the fact 
that the argument of the argcosh in the expression of ip'{m, Sc,ec) given by (j29p is > 1 for all values of m, with two 
point of equalities in m. = and to — TOc . The prescription for the choices of the branches of the arg cosh function 
(that can be changed at the branching point +1) explained after Eq. ([^5]) leaves us with two possible real solutions 
ifi and ip2, which reaches their minimal value in m = and m — rric respectively. These two functions are plotted 
for p = 3 in Fig. [TUl 
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FIG. 10: The two possible eigenstate functions 1^1,2(771, Sc, Cc) at the first order transition point of the p — 3 model. 

This apparent degeneracy of the lowest eigenvalues is however lifted with an exponentially small correction in N. 
Let us denote ap the rate at which this gap closes with N, i.e. 



1 



r(0)/ 



?(0}, 



ap = - hm ^ln( ramjE\"> (s) - E^^" (s)] 



s6[0,l] 

At leading order this rate can be computed from the overlap between the two quasi-eigenvectors 4>i (m) 
and 4>2{'ni) = e"^*^^*^™', as explained at the beginning of this section. We thus have 



. lim -In I (01 102) I 



lim — In 



dm e~^'^i(™'*'='^'=)~^'^2(>n,Sc,ec) 



= mi[ipi{m,Sc,ec) + ip2{m,Sc,ec)] 



(59) 

-Nipi{m) 

(60) 
(61) 



The shape of the sum ipi + ip2 is also displayed in Fig. 1101 It is minimal and constant on the whole interval [0, md; 
indeed the two functions are solutions of Eq. (|29p for the same value of the parameters e, s, and only differ in the 
opposite choice of the branch of the argcosh function for their derivative on [0,r7ic]. As a consequence Op can be 
simply computed by integrating the derivative of </?, i.e. 



a„ = — 
P 2 



dm ach 



Cc + Sc ITT-^ 



(l-Se)VT^^ 



(62) 



This formula, complemented by the values of rric, Sc and Cc as a function of p given in (|14p. is one of the main results 
of the statics part of the paper, giving a very explicit analytical prediction of the exponentially small gap at the 
first-order transition. 

The numerical values of ap thus obtained are displayed in Table HI along with a comparison with the data reported 
by Jorg et al in [3J]. The authors of this paper obtained ap both by exact diagonalization of the matrices iJ^"^ for 
finite N (with an extrapolation in the limit N ^ 00) and by a semi-classical instantonic computation. Our results 
agrees very well with theirs. One can set up an asymptotic expansion of ap at large p, that results in 



In 2 



ar, 



Up 



O 



(63) 



as stated in the Table. The details of this computation are deferred to Appendix [21 The interpretation of ap for even 
values of p shall be discussed in the next section. 



2. The exponentially small gap between the two ferromagnetic phases for even p 



For even values of p the classical part of the Hamiltonian, H{s = 1), is invariant under the reversal of the longitudinal 
magnetization, its groundstate is thus doubly degenerate, with eigenstates fully polarized along the ±z direction. As 
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p 


Tc 


Sc 


m.c 


i^ (diagonalization) [34] 


^ (instanton) [34] 


i^ from Eq. ([62)1 


3 


1.2991 


0.4350 


0.8660 


0.126(3) 


0.1251 


0.1252 


4 


1.1852 


0.4576 


0.9428 


- 


- 


0.2127 


5 


1.1347 


0.4685 


0.9682 


0.270(3) 


0.2686 


0.2680 


6 


1.1059 


0.4749 


0.9798 


- 


- 


0.3057 


7 


1.0873 


0.4791 


0.9860 


0.335(3) 


0.3335 


0.3329 


8 


1.0743 


0.4821 


0.9897 


- 


- 


0.3535 


9 


1.0647 


0.4843 


0.9922 


0.370(3) 


0.3699 


0.3695 


13 


1.0426 


0.4896 


0.9965 


0.410(3) 


0.4105 


0.4093 


17 


1.0318 


0.4922 


0.9980 


0.431(3) 


0.4315 


0.4306 


21 


1.0253 


0.4937 


0.9987 


0.445(3) 


0.4445 


0.4437 


31 


1.0168 


0.4958 


0.9994 


0.462(3) 


0.4623 


0.4618 


p ^>- oo 


1+i 


1 1 

2 8p 


1-^ 


1 1.15 

2 p 


- 


1 ^2 1 

2 12 loK 2 p 



TABLE I: Exponential rate of decay of the gap between the groundstate and the first excited state at the first-order transition 
for the p-spin model, divided by In 2 to ease the comparison with the results of [3J]. The last column is our result, computed 
from Eq. (|62p . The fifth column is the extrapolation from finite A'^ exact diagonalization 34], and the sixth one results from 
an instantonic computation 3^. Thermodynamic parameters Fc = (1 — Sc)/sc, Sc and rric of the system at the critical point 
are also given. The last line gives an equivalent of these quantities at the leading order in 1/p in the large p limit. 



soon as the transverse field is switched on, i.e. for s < 1, this strict degeneracy is lifted. However in the ferromagnetic 
phase, i.e. for s > Sc, this lifting is weak, and the gap between the groundstate and the first excited state is 
exponentially small in iV, of the form e"^^''^'*-' at the leading order. We shall now compute this rate (3p{s), following 
essentially the same lines as in Sec. IIIIG II A similar study for p = 2 can be found in |53 |. 

In the ferromagnetic phase of even p models the infimum in the definition (jlOl) of the groundstate energy is reached 
in ±TO*(s), where the spontaneous longitudinal magnetization to*(s) is solution of Eq. ([TT|) . Hence the argument of 
argcosh in Eq. (j29p is > 1 for all values of m, touching 1 in ±171^.(3). One can thus construct two solutions ip± of 
Eq. (P^ . that vanish in ±m^,{s). An example for p = 2 is displayed in Fig. [TT] As above one obtains the rate f3p{s) 
by computing the overlap between these two quasi-eigenstates. This yields 



/3p(s) = 2 



m* (s) 



d?7iach 



— 7n^. (s) 



.is) 



s mr 



(1 — s)\/l — 1Tl^ 



This formula compares very well with the results of exact diagonalization for p 
of Fig. HH 



(64) 



2 ^56j |. as shown in the right panel 
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FIG. 11: Left panel: the two eigenstate functions ip± for s = 0.4 and p = 2, at the groundstate energy, and their sum. Right 
panel: exponential scaling of the gap /32(s) between the two ferromagnetic solutions for p = 2. The solid black curve has been 
obtained from Eq. (|64p . the symbols are the results of exact diagonalization extrapolated in the limit N ^ oo [5g], the red 
dashed curve is the leading term in the s ^ 1 limit, see Eq. (|67|l . 
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In the classical limit s — )■ 1 the gap vanishes for all values of A'', in consequence the rate /3p(s) diverges in this 
limit. One can study this asymptotic behavior more precisely. One has to*(s) — >■ 1 and egs(s) — > —1 in this limit, and 
the factor 1/(1 — s) makes the argument of the ach function in (j64p diverge for all values of m inside the domain of 
integration. One can then use the asymptotic expansion ach(y) ^ ln{2y) + o(l) for y large to obtain 

1 .1 



/?p(5) - 2 J ^ dm In ^^—-^== ] ^ - ln(l - s) + /3p , (65) 



where the constant /3p can be expressed in terms of the harmonic number function H{x) = L VA-dt, as 



/3p= l\mln('^-^^=^) =\n2- H (^-) +Ih (^-] . (66) 



'0 

In particular for p = 2 we obtain 



/32(s)---ln(l-s)-l + 21n2 , (67) 



which is also plotted for comparison in the right panel of Fig. 1111 

Another interesting limit case concerns the behaviour of /32(s) around the threshold Sc of the second-order transition 
of the p = 2 model. The rate of exponentially small splitting has to vanish in this limit, since the groundstate of the 
paramagnetic phase is no longer quasi-degenerate. More precisely, using the asymptotic behaviors m^{s) ~ 3^/s — Sc 
when s ^ s^ and ach(l + y) ^ y/2y when y — t- 0+, one can expand the expression ([M|) of /32(s) and obtain after a 
short computation that /32(s) = 9(s — Sc)'^^^ + 0((s — Sc)^^^)- This exponent 3/2 was first predicted in [30] on the basis 
of an adaptation of Finite Size Sca ling to mean-field systems (and found also in [57| from the scaling of singular finite 
N corrections). The argument of [30| leads indeed to the value fmfrfc, where z^mf and dc are the mean-field value of 
the exponent controlling the divergence of the correlation length and the upper critical dimension of the universality 
class to which the studied model belongs. In the present case v^i = 1/2, as in the 0** theory, but one has to take 
dc = 3: classical models in this universality class have an upper critical dimension of 4, however in the Suzuki- Trotter 
formulation a d-dimensional quantum model is mapped onto a classical model with an additional imaginary time 
dimension of length /3, and thus correspond to a d -I- 1-dimensional classical model in the zero temperature limit. 
From this value of the exponent and the behavior of the gap in the paramagnetic phase (see Eq. ((52])) the authors 
of |30| could deduce the scaling with N of the gap in the critical regime s ss Sc. Let us reproduce here their argument. 
Suppose that the gap -Ei(s) — Eq{s) satisfies a scaling assumption in the double limit A'^ — >■ oo, s — )■ Sc, i.e. 

Ei{s) ~ Eo{s) ^ N--T{{s - Sc)N-') , (68) 

with T a scaling function and x, x' two exponents to be determined. This assumption can agree with the study of 
the ferromagnetic phase only if x' = 2/3, with J-{z) w exp[— 9^^^^] as z — )■ -l-oo. On the other hand, approaching the 
transition from the paramagnetic phase leads to a closing of the (finite) gap as a square root (see Eq. (|52l) '). hence 
T{z) ~ 2\/2(— 2)^^^ as z — > —00 and x = x' /2 = 1/3. The scaling ass ump tion and the behavior of the gap as N~^^^ 
in the critical regime of the p — 2 model were checked numerically in (30[ . 

Let us finally discuss the structure of the gaps between the lowest states of a model with p > 4 even, in the 
neighborhood of its first-order transition. For the choice of parameters (s, e) = (sc, ec), the argument of the argcosh 
function in Eq. (|29p reaches the branching point 1 in — TOc,0 and +mc, one can thus construct three distinct quasi- 
eigenvectors with rate ip{m) vanishing for these three magnetizations. One could think that the reasoning presented 
at the beginning of this section, that reduces to the diagonalization of a two by two matrix, is invalidated. This is 
however not the case, as is best understood by looking at the three lowest levels of the p = 4 model plotted on Fig. [12] 
On the ferromagnetic side of the transition the groundstate (resp. the first excited state) is the symmetric (resp. 
antisymmetric) combination of the ferromagnetic quasi-eigenvectors concentrated on irric, with an exponentially small 
splitting of order exp[—Nf3p{sc)], while the second excited state is the metastable continuation of the paramagnetic 
groundstate. The avoided crossing of order exp[— iVctp] thus occurs between the groundstate and the second excited 
state; as ap — ^/3p{sc), this gap is much more opened than the one between the ferromagnetic states. The fact that the 
first excited state has no level repulsion at the avoided crossing is easily understood from the additional symmetry of 
even p models discussed at the end of Sec. IIII Al the paramagnetic and symmetric combination of ferromagnetic states 
belongs to the sector invariant with respect to the reversal of the longitudinal magnetization, while the antisymmetric 
combination is in the other, disconnected, sector. Note that on the paramagnetic side of the transition the splitting 
of order exp[—N[3p{sc)] occurs between the first and second excited states. 
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FIG. 12: The energy of the three lowest levels of the p = 4 model in the neighborhood of the first-order transition, obtained 
from numerical diagonalization with A'^ = 40. 



3. Exponentially small gaps between excited states 



In the two previous cases we computed the exponentially small splitting between two quasi-degenerate groundstates, 
the ferromagnetic and paramagnetic ones at the first-order transition in Sec. IIII G II and the two ferromagnetic ones 
in Sec. IIII G 21 It should however be clear (see for instance the left panel of Fig. [7]) that exponentially small gaps 
occur not only between the two lowest eigenstates, but also between excited ones. In this subsection we shall explain 
how to adapt the computation in that case, and in the next one we will in particular obtain the exponential rate of 
closing of the gaps encountered by the metastable continuation of the ferromagnetic and paramagnetic phases, which 
shall be a crucial ingredient for the analysis of the annealing dynamics in Sec. IIVI 

In terms of the leading order eigenvalue equation (|29| , an exponentially small splitting between two eigenstates shows 
up as the existence of two distinct solutions fi.2i'm, s, e) of Eq. (P^ that both fulfill the condition mi[^ip{m, s, e)] = 0. 

m 

We shall call 7(5, e) the rate at which this gap closes, i.e. it is at the leading order of the form e"^'''^*''^-'. As previously 
explained this rate is obtained from the scalar product between the two quasi-eigenvectors. The two functions (^1^2 
only differ by a choice of branch of the argcosh function on an interval [mi(s, e), 7712(5, e)]. This implies that their 
imaginary part is the same for all m (see Eq. ((30)) ). hence the scalar product between the two eigenvectors depends 
only on the real part of (/'i,2- This leads to: 



7(s,e) = - lim — In |(0i|(/.2)| = - lim —In 

W-s-oo 1\ N-¥oo I\ 



dm e~^'^i(™''*'^)~^'''2(™''*'^) 



1 rm2{s,e) 

mi^[ipi{m,s,e) + ip2{iTi,s,e)] — - / dmach 

™ 2 J„ii{s,e) 



{l's)VT 



(69) 



Let us now describe the regions in the (s, e) plane where exponentially small gaps occur. The discussion above 
shows that their occurence can be traced back to the number of times the argument of the arg cosh function in Eq. (P^ 
reaches the branching points ±1, in other words the number of solutions m £ [—1, 1] of the equations 



T.P-il~s)^/T 



or e 



-smP + (l-s)VT 



(70) 



A moment of thought reveals that for any value of (e, s) in the authorized range of eigenvalues (defined in ([5T|) ) 
this number is either 2, 4 or 6 (counting twice the marginal case of a branching point touched quadratically and 
not crossed). The first case corresponds to a non-degenerate eigenstate, the two others to exponentially small gaps 
between eigenstates. The frontiers between these domains correspond to the disappearance of some solutions of the 
equations ([70]) . which define implicitly ?ti as a function of s, e. Their boundary can thus be obtained as the limits of 
validity of the implicit function theorem. After a short computation one realizes that they are given by curves e(s) 
of the form ((7(I| . with ni replaced by one of the solutions (stable or instable) of Eq. ([TT|). Let us be more concrete by 
distinguishing between various cases: 

• for p = 2, when s < Sc the spectrum is made of non-degenerate eigenvalues with e G [—(1 — s), (1 — s)]. When 
s > Sc the low-energy part of the spectrum (e € [egs(s), —(1 — s)]) has doubly degenerate eigenvalues with an 
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exponentially small gap between them, while the high-energy spectrum (e €! [—(1 — s), 1 — s]) is non-degenerate. 
These two regimes are depicted in the left panel of Fig. [T31 and agree with the qualitative features described in 
Sec. unci on the basis of the numerical diagonalization (cf. the left panel of Fig. [5]) . For clarity a zoom of the 
latter is presented on the right panel of Fig. [IJl in the neighborhood of the line e = —(1 — s), which is indeed 
the point where the energy splitting of excited ferromagnetic states is no more exponentially small. 



e 





FIG. 13: Left panel: the two regimes in the (s, e) plane for the p = 2 model. Right panel: a zoom of the spectrum obtained by 
numerical diagonalization for A'' = 80 around the line e — —(1 — s) where the gaps between excited ferromagnetic states are no 
longer exponentially small. 

• for p > 3 odd, the spectrum is symmetric under e -^ — e (as explained at the end of Sec. IIII A|) . we shall thus 
describe only its part with negative e. The equation (llip has only tti = as a solution for s < Sgp, with an 
associated energy epm(s) = —(1 — s), while for s > Sgp there are three solutions < mi(s) < to*(s) with energies 
epm(s), ei(s) and efm(s). These three energy curves are drawn on the plots of Fig. [T3]for p — S; the energies 
corresponding to doubly-degenerate eigenstates with exponentially small gaps between them are in the range 
[max[efi„(s), epm{s)], ei(s)] for s > Sgp. On the right panel we have superimposed the spectrum for N = 320, one 
sees indeed that the avoided crossings occur precisely in this regime. All other authorized values of the energy 
correspond to non-degenerate eigenvalues. 

• for p > 4 even, the phenomenology is mixed between the one of the p — 2 and the p > 3 odd cases. Three 
zones are to be distinguished in the (e, s) plane (see Fig. 1151) . One corresponds to the ferromagnetic phase, with 
doubly quasi-degenerate eigenstates of opposite magnetizations, for s > Sc and e € [egs(s), epni(s)]. In the area 
s > Ssp, e € [max[efin(s), epiii(s)], ei(s)] there are three valid solutions of Eq. (j29p . As explained at the end of 
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FIG. 14: Left: the three areas in the negative energy part of the spectrum of p = 3. Right: a blow up of the data from 
numerical diagonalization in the region of exponentially small gaps for A^ = 320. 
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FIG. 15: Left: the three areas in the spectrum of the p = 4 model. Right: avoided crossings between paramagnetic excited 
states and quasi-degenerate ferromagnetic excited states, for p — 4 and A'^ = 40. 

Sec. IIII G^ avoided crossings in this area, that are of order e"^'''*^''''^', only occur between continuation of levels 
coming from the paramagnetic zone and combination of ferromagnetic quasi-eigenvectors that have the same 
parity under the reversal of the longitudinal magnetization. The splitting between symmetric and antisymmetric 
combinations of ferromagnetic states is much smaller, of order e~'^^'^^^'^' . This phenomenon comes from the 
additional symmetry of even p models discussed at the end of Sec. IIII A[ and is illustrated on the right panel of 
Fig. [ini The other authorized regime in the (e, s) plane leads to single solutions of Eq. ([2^ . 



4- Exponentially small gaps encountered by the metastahle states 



In view of the application of these computations to the annealing dynamics in Sec. IIVI the most important case to 
consider is the size of the gaps encountered along the metastable continuations of the paramagnetic and ferromagnetic 
groundstates. We shall thus define 7pm(s) = 7(5, epm(s)) for s e [sc, 1] and 7tm(s) — 7(5, efin(s)) for s e [sgp, Sc]. These 
quantities are plotted for p = 3 in Fig. 1161 along with an example of the functions if 1,2 involved in the computation 
of 7pni for one value of s. The numerical evaluation of these quantities is easy thanks to the explicit expression (|69p . 
One can also perform analytically some expansions around special values: 

• 7fin(s) vanishes at Sgp as 7p(s — Ssp)^'**, with the prefactor expressed as 



Ip 



6%/2 (p-1) 






l+P 



(p-2) 



b-1) 



(71) 



The exponent 5/4 is in agreement with the reasoning of [301 recalled in Sec. IIII CTS] Indeed the spinodal transition 
is in the universality class of cubic field theories, with the upper critical dimension (taking into account the 
imaginary time direction) dc = 5, and the mean- field value of the critical exponent for the divergence of the 
correlation length i^mf — j . One can also adapt the Finite Size Scaling argument of [3^ to predict that for large 
but finite values of N the gaps encountered in the neighborhood of the spinodal should scale as 7V~*/^^. This 
follows from a scaling hypothesis of gaps of the form given in Eq. (j68l) . combined with the exponent x' = 4/5 
obtained above from the limit s — >■ s+ , and the closing of the finite gaps in the limit s — >■ Sgp, argued to occur 



with an exponent 1/5 at the end of Sec. IIII Fl 

On the other hand the vanishing of 7pm when s — >■ 1 is non- universal (i.e. depends on p), one finds indeed 



7pm(s = 1 -(5) ~ 7p(5p-2 



7p 



— p— / dx xyl — xP^'^ 
2F^ Jo 



(72) 



• In the neighborhood of Sc the behavior of 7fm and 7pni exhibit a singularity of the form (s — Sc) ln(s — Sc), more 
precisely 



7pm (sc + 6) -^ ap + r]pSln{6) , 7fm(sc - 6) -^ ap + fjpSln{6) 



(73) 
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FIG. 16: Left panel: the exponential rates 7pm(s) and 7fm(s) of the gaps encountered along the metastable continuation of the 
paramagnetic and ferromagnetic groundstates, for p = 3. Right panel: the functions (^1,2 yielding 7pm(s = 0.45). 



where the constants 'qp and rjp are given by 



Ip 



(p-l)p-i 
p"^{p— 2)t~ 






p-i 



rip = V(p-l)(p-2)?7p 



(74) 



IV. QUANTUM ANNEALING OF THE MODELS 



This section is devoted to the study of the annealing of the models whose static properties were considered above, 
and is organized as follows. We shall first (in Sec. lIV Al) define precisely the dynamics and the quantities of interest to 
be studied in this context. Then we will review the phenomenology of the simple, two-level, Landau-Zener problem 
in Sec. IIVBI and discuss on this basis the expected features of the annealing dynamics. A further simplified model 
is introduced as an aside in Sec. llVCi that will be used as a benchmark for the comparison with numerical results. 
The actual computations and results will then be presented in two sections, divided according to the scaling of the 
annealing time with the system size; in Sec. lIVDJ we shall consider annealings on exponentially large times, while in 
Sec. I IV El we will study the behavior of the dynamics when the thermodynamic limit is taken for a finite annealing 
rate. The main results of Sec. IIVDI and IIVEI are presented for odd values of p > 3 for which the phase-transition is 
first order and not mixed up with the quasi-degeneracy of the ferromagnetic states. We briefly comment in Sec. IIVFI 
on the behavior for even values of p, in particular p = 2 (the Curie- Weiss model) whose annealing was already studied 
in MM. 



A. Definitions 

As sketched in the introduction the quantum annealing procedure, or quantum adiabatic algorithm, aims at finding 
the groundstate of some final Hamiltonian H{ via an interpolation from an initial Hamiltonian Hi whose groundstate 
is easy to construct. The system evolves from time i = to i = T, the total running time of the algorithm, according 
to the Schrodinger equation with an Hamiltonian interpolating (for instance linearly) between Hi and iff. In terms 
of the reduced time s — t/T e [0, 1], this reads 



Td^ 



-\Ms)) ^ H{s)\Ms)) 



H{s) ^{l- s)Hi + sH{ 



(75) 



with the initial condition that 10^(0)) is the (normalized) groundstate of Hi (we set h — 1 from now on). The outcome 
of the algorithm for an annealing time T is thus the final state \(J)t{1)), which ideally, if T is much larger than the 
adiabatic time, is close to the groundstate of Hf. 
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It is however interesting, in particular for the approximabihty issues mentioned in the introduction, to study this 
procedure also for T smaller than the adiabatic time. We shall quantify the deviation from adiabaticity by computing 
the final energy density defined as 

eUT,N) = l(0^(l)|i/f|0j,(l)) , (76) 

and comparing it to the groundstate energy density Cgg of the final Hamiltonian iJf : the residual energy density is thus 
Bros = Cfin — Ggg. Another relevant energy density to compare Cfln to is the trivial one achieved when the interpolation 
time vanishes, i.e. when one computes the average energy of the final Hamiltonian with respect to the groundstate 
of the initial one: etriv — ;^(0T(O)|i?f |</>t(0)). Indeed egain — etriv — efin is the gain in energy density that is achieved 
by the evolution during the time T . Note that with the normalization we chose for the models one has Cgs = — 1 and 
etriv = 0. 

Our analytical results will all be obtained in the thermodynamic limit iV — > oo, but with two different scaling of T 
with N that shall be distinguished typographically. If T is kept fixed when N diverges we shall denote 

efi„(T)- lim efi„(r,iV) ; (77) 

this regime will be studied in Sec. IIVEI On the other hand if T scales exponentially with N (as in Sec. IIVPP we call 
T this exponential rate and define 

efin(r)= lim efi„(^ = e^^Af) . (78) 

A*— >oo 

We shall argue in the following that, as far as intensive quantities like the energy density are concerned, these two 
regimes are the only relevant ones for p > 3 (see Sec. lIVFl for a discussion of the different case p = 2), i.e. polynomial 
scalings of T with N are just limiting cases of the two regimes above. Note also that in the thermodynamic limit, for 
both regimes, the quantum fiuctuations of the final energy density are neglectible, hence a description in terms of the 
average energy only is meaningful. 

The definitions above are valid for any choice of the inital and final Hamiltonians. From the point of view of 
potential applications they are of course most interesting when Hf has a groundstate that is a priori hard to find and 
when it is easy to prepare the system in the groundstate of Hi. In the following we shall consider the dynamics of the 
annealing of the models whose statics were studied in the first part, that is use a p-spin interaction and a transverse 
field as initial and final Hamiltonian. Obviously neither of these Hamiltonians has a groundstate which is hard to find, 
hence they can only be considered as toy models for the application of the quantum adiabatic algorithm. However 
they share some properties (first-order transitions, metastability, spinodals) with more realistic random combinatorial 
optimization problems (iTl |l8| , while being much easier to study both analytically and numerically. Because of this 
unrealistic character both choices of the transverse field as Hi and the ferromagnetic interaction as H{ or viceversa 
are equally relevant, and it will be very instructive to consider these two types of evolution. Let us define them more 
precisely: 

• The annealing towards the ferromagnet corresponds to the choice Hi = —Nrfv^, Hi = —N(fh^Y , i.e. 
H{s) = — (1 — s)Nfh^ — sN{fh^)P is precisely the Hamiltonian ([2]) studied in the first part of the paper. 

• The annealing towards the paramagnet corresponds to the reverse choice Hi = ~N{fh^Y , Hi = ~Nfh^ , in other 
words the evolution with the Hamiltonian ^ is made with s decreasing from 1 to 0. To avoid confusion we 
shall denote u = 1 — s instead of s the reduced time in this case, i.e. study the following equation: 

^^\Mu)) = [-iV(l - u){mY - Num^UTiu)) . (79) 

Note that in all these cases the groundstate of the initial Hamiltonian Hi belongs to the fully symmetric sector of 
maximal spin. The instantaneous Hamiltonian H{s) is block diagonal with respect to the spin decomposition for all 
values of s, in consequence the state \4>t{s)) remains in the maximal spin sector i^ = all along the evolution. 

B. Finite duration Landau-Zener problem and its expected consequences 

The Landau-Zener problem |58l . |59| is the simplest example of a quantum evolution with a time-evolving Hamilto- 
nian. It involves two levels of linearly varying energy with a fixed coupling between them: 

(80) 
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The initial condition is given by ipi{t — >■ — oo) = 1, i.e. the system is initially in its groundstate. The probability of 
transition to the excited state after an infinite time can be computed exactly (see |60l - l62| for modern derivations) and 
yields Poxc = lim 1-01(^)1^ = exp(— Trae^). It is thus a function of the product between the square of the minimal gap 

i— f+C30 

e at i = and the velocity a of variation of the energies of the levels. 

Variations of the Landau-Zcner model that account for a finite duration of the interaction have been studied in 
great details in 63', 643 • Consider for instance an evolution with a reduced time s S [0, 1] corresponding to a total 
physical time T, with two levels that have an avoided crossing at s = 1/2: 




—a{s — 




(81) 



The probability Poxc(a, e,?") that the evolution starting from the groundstate at s = leads to the excited state at 
s = 1 can be expressed in terms of special functions [63| and simplified in various asymptotic limits according to the 
relative ordering of a, e and 1/T. In the present context the relevant regime corresponds to a fixed, e — ^ and T — >■ oo. 
Then Pckc has a scaling form if T diverges as e^^, more precisely 

lim Pcxc(a, e,T — ae^^) — exp[— vra/a] . (82) 

If one further assumes that both T and e scales exponentially with a large parameter N, according to e{N) — e~'^^ 
and T = e^^ , then the probability of excitation reduces to 

lim Fexc = ^(27 - t) , (83) 

with 6{x) the Heaviside step function, i.e. on this scale either the evolution is sufficiently slow and the system follows 
adiabatically the groundstate or it is too fast and with probability 1 the system goes into the excited state. 

Let us now explain the intuitive picture for the dynamics of the p-spin ferromagnetic model (with an odd value 
of p > 3) in the large N limit, that arises from the combination of the study of this two-level problem and of the 
results of Sec. IIIII (a similar reasoning can be found for instance in [65|)- Consider first the annealing towards the 
ferromagnet, for a large evolution time T, starting from the groundstate at s = 0. As long as s < s^ the gap between 
the groundstate and the first excited state remains finite, hence for times sufficiently large (but independent of the 
system size), it is expected that the system will remain in the instantaneous ground state. There occurs at Sc an 
avoided crossing with an exponentially small gap of order e"^"". Transposing the results of the two-level problem, 
two cases have to be distinguished. If the evolution time is exponentially large, T = e^^ , and if r > 2ap, then the 
system follows adiabatically the groundstate at the avoided crossing, and continues on the instantaneous groundstate. 
Otherwise the system is in the first excited state just after the crossing, i.e. on the metastable continuation of the 
paramagnetic groundstate. We have seen that this state encounters a series of avoided crossings, that lead to gaps 
of order e~^'*'p"*^*''. Let us assume that all these avoided crossings are independent, and can be treated as in a 
two- level problem. Then, if the evolution time is T = e , one is led to the conclusion that the system will remain 
in the metastable groundstate until the value stum such that r = 27pni(sturn), and from thereon follows the excited 
ferromagnetic state that crossed the paramagnetic metastable state in Sturn- As there is no spinodal limit for the 
metastable paramagnet 7pm (s) > for all s < 1. Hence for an evolution on sub-exponential times T the system 
should follow the paramagnet until s = 1, which leads to a vanishing energy density gain with respect to the trivial 
one. 

A similar reasoning in the case of the annealing towards the paramagnet reveals a richer phenomenology. For an 
exponentially large annealing time e^"^ with t > 2ap the groundstate is followed during the whole evolution. If 
T < 2ap the metastable ferromagnetic state will be followed until the turning point uturn where r = 27fni(l — wturn), 
then the system follows the paramagnetic excited state that rejoins the metastable ferromagnet at the turning point. 
There is however an important difference with respect to the reverse direction of annealing: here the ferromagnet has 
a spinodal limit of metastability. Hence an evolution on an exponentially long time e^"^, but for arbitrarily small 
values of r, yields a non-trivial (negative) energy density. In addition the regime of large but sub-exponential T can 
be expected to be much richer than in the previous case: the ferromagnet will be followed for u < Ugp = 1 — Sgp, but 
for subsequent times this analysis in terms of level crossings can give no clue. 

In this reasoning we have assumed that the various level crossings can be treated independently one from the others, 
and apply to each of them the results of a simple two-level problem. Arguments in favor of this assumption can be 
found from static [S^l and dynamical [62, |6J| considerations: as can be seen on the drawings of the spectrum (see 
for instance the right panel of Fig. 15]), an avoided crossing with an exponentially small gap affects notably the two 
colliding levels on an interval of s which is also exponentially small. On the other hand two successive crossings are 
located at values of s which are distant of order 1/iV. Similarly in the dynamical case the "duration" of a crossing 
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(as defined in ^Jl) should go like exp((— 7 + t)N), and therefore the influence of a crossing should spread on a range 
of s of order at most exp(— 7A^). 

In the following sections we shall present the explicit results obtained from this reasoning, and compare them with 
the results of numerical integration of Schrodinger's equation for finite values of N. In particular we will test the 
assumption of independence of the different crossings. The numerical results for the annealing of the p-spin model 
having strong finite-size corrections, we shall first introduce a simplified model that shares some of the properties of 
the p-spin model but with smaller finite-size effects. 

C. A further simplified model (the p — >^ 00 limit) 

We shall introduce here a simplified version of the models under study, first defining it formally and discussing 
afterwards its relationship with the main models of the article and with previous works. 

We consider an interpolating Hamiltonian Ht{s) acting on the fully symmetric subspace of dimension A^ + 1. It 
is given by -ff»(s) = (1 — s) J — Nsfiv^ ^ with the initial Hamiltonian H^ = J defined by its matrix elements in the 
x-diagonal basis. 



{m;0\J\m';0), = -ND„^D„ 



J-^rn 



1 



N 



2N\j^l±rn 



(84) 



It is thus a matrix of rank 1, with a single eigenvalue equal to — A^ and all other eigenvalues equal to 0. The spectrum 
of H,{s) is presented on the left panel of Fig.[T71 As J is of rank one the spectrum of Ht{s) is essentially equal to the 
one of —Nsfh^ (see for instance [66] for general results on low rank perturbation theory). There is however a major 
difference: the isolated eigenvalue of energy density e = — (1 — s) is continued as a metastable state for s E [1/2, 1], 
with exponentially small avoided crossings of order e^^'''*^'*). The computation of 7,(s) is presented in Appendix IB 1[ 
and yields the explicit formula 



7.(s) = -^lns + ^^ln(2.s-l); 



(85) 



this function is plotted on the right panel of Fig. 1171 For the reasons explained above one expects that the annealing 
of the model on sub-exponential time scales yields a vanishing final energy density, as the metastable continuation of 
the groundstate exists until s = 1; this is confirmed by the analysis presented in Appendix IB 21 Exponentially slow 
annealings with < r < ln2 (i.e. T ^ 2^) should however reach a non-trivial negative energy density, larger than 
the one of the groundstate but smaller than the trivial one, egg < efin(T) < etriv in the notations of Sec. IIV A1 




FIG. 17: Left: spectrum of the operator H,{s) as a function of s, for TV = 80. In the thermodynamic limit the metastable 
continuation of the groundstate at s < 1/2 exists until s = 1. Right: The exponential rate of closing of the gaps at the avoided 
crossings 7. (s), defined in Eq. (|85p . 



This simplified model is actually (almost) the p — )■ cx) limit (with p odd, and the limit on p taken before the limit 
on N) of the models studied in the main part of this paper (and was discussed in these terms in [3J|): in the z-basis 
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J is diagonal, with matrix elements z{m; 0\J\m; 0)^ — —NSm.i, to be compared with lim z{'m', 0| — N{fh'^)P\m; 0)^ — 

—N5m,i + NSrri,-! if the limit is taken with p odd. If one focuses on the low-energy part of the spectrum, one can 
view the evolution of the simplied model as the evolution towards the paramagnet of the p-spin model in the large 
odd p limit. 

Another justification for the introduction of this simplified model can be given as follows. Assume that one is given 
an arbitrary Hamiltonian H{, diagonal in the computational basis of the 2^ classical configurations of spins, as an 
optimization problem, and that the problem is to be solved without using any information about the local structure 
of these energies in the configuration space. Then the most natural starting Hamiltonian Hi for an interpolation is 
the one connecting any two configurations of the Hilbert space with equal probability, 

i^i = ./=-^El^)(^'l = -^l^)(^l' ^ith |X) = ^E|a), (86) 

where the normalization chosen is such that J has one eigenvector \X) with eigenvalue —A'' and 2^ — 1 eigenvectors 
with eigenvalue 0. Let us denote {i?Q}ag[i.M] the distinct energies of Hf, da the number of configurations a on which 
Hi takes the value £"„, and H the Af -dimensional Hilbert space generated by the symmetric combinations of the 
states of a given energy: 

H = span{|a)}, 1") = -^ E 1^) ' (8^) 

Then the ground state \X) of J is in H, and so is, for any s, the vector \(J)t{s)) obtained by the evolution according 
to the Schrodinger equation with H{s) = (1 — s)J + sH{ as interpolating Hamiltonian. The dynamics can thus be 
studied in the symmetric subspace Ji, in which the matrix elements of H{s) are given by 



{a\H{sm = sSa.pEa - (1 - s)N^^^ . 

The simplified model defined at the beginning of this section is thus a representative example of this more general 
construction, in which we chose M = N + 1, with equally spaced levels Ea between — A^ and -l-A^, each with a binomial 
degeneracy. The quantum annealing with such an unstructured Hamiltonian J has been studied in [67] . In the context 
of Grover [631 search problem (i.e. with a golf course potential H{ having a single low energy level), it was shown 
in [69| that a modification of the annealing procedure could reproduce Grover's quadratic speedup. By slowing down 
the interpolation in the neighborhood of the avoided crossing one can indeed reduce the adiabatic time to 0{2^^'^). 

D. Annealing on exponentially large times 

1. The simplified model 

Let us compute the final energy efin(''') after an exponentially long annealing of duration T = e^'^, for the simplified 
model of Sec. llVCi following the reasoning of Sec. IIVBI The turning point Sturn up to which the metastable state is 
followed is given implicitly by 27,(sturn(''')) — t, where the expression of 7, is given in Eq. ([55)1 (if r > 27,(1/2) — In 2 
we set Sturn(''') = 1/2)- The final energy at the end of the annealing is then given by the continuation of the state that 
crosses the metastable state at Stum- For this simple model where energy levels are at leading order linear functions 
of s except at the crossings, this yields: 

efin(r) ^ - ''-'7^ = 1 - -rf-- , (89) 

Sturn (t) 7. (t/2) 

where 7"^ is the functional inverse of 7,, with the convention that j^^{z) = 1/2 if z > i^. 

A comparison of this analytical prediction with the results obtained by numerical integration of the Schrodinger 
equation (see Appendix[C]for details on the procedure we used) is presented in Fig.[T8l The left panel displays the result 
of Eq. (|M|) along with curves efin{T — e^'^, A^) obtained numerically for some finite values of Af. We extrapolated these 
results in the N ^ 00 limit with finite size corrections of the form efin{T = e^'^ ,N) = a{T) +b{T)^-^j^ + c{t) j^^ + o{l / N) , 
a form that can be expected to arise because of polynomial corrections to the exponentially small gaps; the inset shows 
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the very good quality of such a fit already for small values of N. The extrapolated curve a{T) agrees with the analytical 
prediction 1^ within 1%. 

The right panel of Fig. [T8]provides a further confirmation of the analysis in terms of independent two level Landau- 
Zener problems. The black symbols with error bars represent the quantum average and standard deviation of the 
instantaneous energy, 



eris) = l(^j,(s)|i?(s)|0T(s)) , 



aris) = l^^{Ms)\H{sr\M^)) - (^{Mmis)\Ms)) 



(90) 



computed numerically during an evolution with N = 64, T — e^"^ for r — 0.2. One observes indeed that the average 
instantaneous energy follows the metastable groundstate across several crossings, until the turning point after which 
it follows adiabatically the levels crossed there. The standard deviation is almost constant in time except around the 
turning point where it grows slightly, reflecting the fact that for finite N a few levels (those with gaps close to T~-^/^) 
get populated. The independence of the crossings is even more apparent in the inset, which shows that the slope of 
eris) jumps significantly for three values of s that correspond precisely to the locations of avoided crossings. 
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FIG. 18: Left panel: final energy density for the evolution of the simplified model as a function of r = {logT)/N. The sofid fine 
is the analytic prediction (|89[) . The symbols are the results of the integration of Schrodinger equation for A^ — 16, 36, 64, 128, 
and an extrapolation to N —^ oo using corrections in {In N)/N and 1/N. The insets shows this fit for r = 0.1, with a fitting 
function of the form f{N) = a + bln{N)/N + c/N, with the results of a best fit on TV > 20 given by a = -0.438, b = -1.463, 
c = 2.926. Right panel: black symbols and error bars represent the average and standard deviation of the instantaneous energy, 
for the evolution of the simplified model with A'' = 64, r = 0.2. The red lines correspond to the spectrum of H,{s). The inset 
is a zoom around the turning point, the arrow on the right is the prediction of Eq. (I89|l for the final energy density in the 
thermodynamic limit. 



For completeness let us state the asymptotic expansions of enniT) around t — and r = In 2, that are easily 
deduced from the behavior of 7,(s) in s = 1 and s = 1/2, respectively, and read 



efin(T) -- -V2t , 



efi„(T = ln2-(5) ^ -1 + 2——- 
s^o+ ln(l/(5) 



(91) 



2. The annealing towards the ferromagnet 



We now follow the same reasoning for the annealing of the p-spin model towards the ferromagnet. The turning point 
Sturn is given by Sturn(''") — 7~i(T/2), where the function 7pni was computed in Sec. IIII G 41 and plotted on the left 
panel of Fig. [1^1 We adopt again the convention that Jpn-^iz) = Sc if 2 > 2ap. As already mentioned in Sec. IIII El the 
computation of efin(T) is then completed by an iso-density argument: for s > Sturn(''') we assume that the evolution 
follows adiabatically the eigenstate that made an avoided crossing with the paramagnetic metastable state at Sturn- 
Because of the absence of any level crossing in this regime the number of eigenvalues below the one whose energy we 
want to follow is constant, by definition. Hence efin is fixed by the condition 



2?o(sturn(T), -(1 - Sturn (t))) = 2?o(l, efi„(r)) = 



1 - (-efi„(r))VP 



(92) 
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where the integrated density of states X'oCsje) is given in Eq. P5|) . and the last equahty follows from its explicit 
expression when s = 1, p is odd and e < 0. This prediction is displayed in the left panel of Fig. 1191 along with results 
of the numerical integration of the Schrodinger equation for finite N. The finite-size effects in these results are much 
larger than for the simplified model. The extrapolation towards iV — >■ cc was done searching the value of r that 
corresponds to a given final energy density instead of the contrary (see the caption of Fig. [TH]for details), and gives 
a satisfactory agreement with the analytic prediction. 
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FIG. 19: The final energy density for the evolution on exponentially large times, as a function of r = {\ogT)/N. The left (resp. 
right) panel corresponds to the annealing towards the ferromagnet (resp. paramagnet) for p — 3. The solid lines at the bottom 
are the analytic prediction from Eq. H92|) for the left panel, and Eq. H96|) for the right one. The other lines with symbols results 
from the integration of Schrodinger equation for various finite sizes. The black symbols on top of the analytic predictions are 
the extrapolation in the N ^ oo limit of the numerical results. The latter was performed by fitting, for various values of e, the 
(exponential) time t{N) such that efln(r = g^'^(^)^ jv) crossed e, with a fitting function of the form t{N) = r + cst^^^ -fcst-^. 
We show in the insets the details of the fit for e = —0.5 (left panel) and e = —0.945 (right panel). 

The final energy density vanishes in the r — > limit, as the paramagnetic metastable state has no spinodal and can 
thus be continued until s = 1. A more precise asymptotic statement can be obtained by studying the behavior of the 
integrated density of states close to s = 1, namely 



1 i i-p 1 f i-p 

2?o(l — S, —S) ^ dp6p , dp — 2 p / dxx p acos (1 — x) , 

2 Svr Jq 



(93) 



for odd values of p. Combining this expansion with the one of 7pm (s = 1 — (5) stated in Eq. (|72p yields 

efi„(r) ^~o~Sp'^^ ' ep:^2^dPj^ . (94) 

On the other hand the behavior of efin(''') for (exponential) times slightly smaller than the adiabatic time r = 2ap 
corresponds to the limit where Stum — ^ s^ ■ One can thus invert the expansion (|73p for the behavior of 7pni around Sc 
to obtain the behavior of Sturn(''' ~> 2ap). Noting that 2?o(s, ^ (f — s)) has a finite derivative with respect to s in Sc, 
one obtains finally after the simplification of various constants: 



efi„(T = 2ap - (5) - -l + 2p— — -— , 
(5->o+ ln(l/o) 



(95) 



whose form is similar to the one found for the simplified model in Eq. (|9T 



3. The annealing towards the paramagnet 



The annealing towards the paramagnet can be treated along the same lines. We recall that in this case the 
interpolation parameter is u — 1 — s. The evolution follows the metastable ferromagnet until the turning point 
Mturn G [f — Sc, f — Ssp] such that T = 27f,„(l — Wturn), with 7fm(s) the function computed in Sec. IIII G 41 and plotted 
on the left panel of Fig. 1161 The isodensity argument for the continuation of the evolution in the regime u > Uturn 
reads then 



2?o(l - '«tum(r), efni(f - Uturn(T))) = 2?o(0, efin(r)) 



f +efin(T) 



(96) 
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the last equality being the consequence of the equidistance of the paramagnetic levels when u = 1. A comparison of 
this analytical prediction with numerical results is shown on the right panel of Fig. 1191 The agreement of the large N 
extrapolation with the prediction of Eq. (|96)) is again satisfactory. 

The small r limit of this regime yields a non-trivial energy density, because the ferromagnetic metastable state has 
a spinodal limit of existence. It is given by the continuation of the paramagnetic state that goes to the spinodal point, 
and from the formula above reads: 



Cfin = lim efi„(r) = -1 + 2Po(ssp, e^ 



(97) 



We report the value of these energy densities for some values of p in Table |TT1 For large p, using the asymptotics 
of (J13p . and the fact that the energy of paramagnetic levels become linear functions of s in this limit, one gets the 
asymptotic behaviour Cfin ^ — 7=, s^ form that agrees very well with the data in Table |TT1 

p->oo vP 

The correction of next order in r is obtained from the asymptotic expansion of 7fni around Sgp, given in Eq. (|71L 

1 - Ss 



which converts into Ut 
d 



M„ 



ds 



urn v"^/ 

2?o(s,ef,n(s)) 



(t/27p)'*/^. Let us define the positive constant 
1 



2tt{1 - Ssp) 



dm ■ 



™?p 



^(1 - s,p)2(l - m2) - (csp + Ssp mPy 



(98) 



where m' is defined as the negative value of m where the square root vanishes (one can notice that Mp = ^(0), where 

Then expanding in Eq. (|96p one obtains the final energy density 



G{z) is the scaling function defined in Eq. (|57] 
behaviour as 



efin(T) 



efin - 2Afp 



2% 



(99) 



The opposite limit of quasi-adiabatic times yields exactly the same formula for the energy density as in the simplified 
model (which is indeed its p — >■ oo limit), i.e. 



efin(r = 2ap - S) 



-1 



(5-i.0+ 



ln{l/6) 



(100) 



p 


Tsp 


Ssp 


rrisp 


efln 


3 


1.5 


0.6 


0.7071 


-0.9302 


4 


1.540 


0.6062 


0.8165 


-0.8259 


5 


1.624 


0.6189 


0.8660 


-0.7861 


7 


1.812 


0.6443 


0.9129 


-0.6881 


9 


1.994 


0.6660 


0.9354 


-0.6187 


13 


2.325 


0.6993 


0.9574 


-0.5256 


21 


2.884 


0.7426 


0.9747 


-0.4211 


31 


3.462 


0.7759 


0.9831 


-0.3500 



TABLE II: Final energy for an annealing of the p-spin model towards the paramagnet, in the limit of "small exponential" times. 
The thermodynamic parameters of the system at the spinodal point are given by (|13|) . 



E. Annealing on constant times 



We now turn to a study of the dynamical properties of the previous models on time scales not growing exponentially 
fast with the size of the system. From the analysis of Sec. lIVBI we expect that for the simplified model and for the 
annealing towards the ferromagnet the final energy density vanishes on such time scales, because the metastable 
branch with exponentially small avoided crossing exists until 5 = 1. This is confirmed for the simplified model by 
a technical analysis that is deferred to A ppe ndix IB 21 The annealing towards the ferromagnet can be treated via a 
semi-classical dynamical analysis |39l l40l |43|. and this will confirm its triviality on constant time-scales. The same 
semi-classical analysis will on the other hand reveal a rich structure for the annealing towards the paramagnet on 
finite time scales. 



32 

1. Semi-classical dynamics for the annealing towards the ferromagnet 

Let us decompose the vector \(J)t{s)) on the z-diagonal basis as 

\Ms))= Y. <^T(m,s)|m;0), . (101) 

The Schrodinger equation (I75p is equivalent to a set of coupled equations for these coefficients, 



i d(j)Tim,s) „ , , , (1-s) / .2, , / 2 

= -8 11%^ (j)T[m,s) - — Wl -m^ + — (1 -to) 0T m + — ,s 



NT ds -r.y , , 2 V N^ ' ^ \ N 



which is the analog of ((27|) in the stationary case. The semi-classical dynamic Ansatz is 4>T{m, s) = Q--^VT{m,s) ^ which 
yields in the large N limit, with T fixed, the evolution equation 



I diprim, s) ^ _g ^p _ (^ _ gW^ _ ^2 cosh (2y'r(m, s)) , 
T OS 



(103) 



where the prime denotes the derivation with respect to m. This corresponds to (j28p with the replacement e —>■ — y-^- 
The initial condition is the groundstate of the pure transverse field, it is thus given by (^T('7i, 0) = ipoirn), with ipo 
defined in Eq. (j34]). This partial differential equation is rather difficult to solve numerically. One can however make 
a further analytical simplification. 

The computation of physical observables that are diagonal in the to^ basis only requires the knowledge of the location 
of the minimum of the real part of the large deviation function ipx , that we shall denote qxis) — arg min„i 5R(pT ("t-, s). 
In particular at the end of the evolution the final energy is given by efin(T') = — (j't(I)p. It turns out, as explained 
in [43| . that it is possible to write a closed system of two differential equations on qris) and its conjugate momentum, 
5t(s) = —dm'^VTiqT{s),s). The evolution equation (I103P implies indeed 

T^lTis) = -^niqT{s),qT{s),s) = 2(1 - .s)v/l - qris)^ sin(2qT(s)) , 



l^qris) = -|-H(gT(s),$T(s),s) = spqrisf-' - (1 - s)^=Mkl==cos{2qT{s)) . 



(104) 
T ds dq ^l -qris)^ 

These are Hamilton equations of classical mechanics, with an Hamiltonian 



niq, q, s) ^-sqP-{l- s) ^/ 1 - q^ C0s{2q} (105) 

obtained from the differential operator on the r.h.s. of (jl03p by the canonical substitution m ^ q, i-^ -^ q- For the 
sake of completness we explain in Appendix |D] the derivation of (jl04p from (jl03p , along the same lines as in [43[ ; note 
also that similar semi-classical equations can be obtained for fermionic models within the time-dependent Gutzwiller 
approximation 70]. 

One can in addition show that the average instantaneous energy of the evolution according to the Schrodinger 
equation is precisely equal to the classical Hamiltonian, namely 

lim ^lcj,T{s)\H{s)\Ms))^H{qT{s),qT{s).s) . (106) 

Let us now conclude on the validity of the analysis of Sec. IIVBI i.e. that for annealing times T that are constant 
in the thermodynamic limit the final energy efin(T') vanishes. The initial condition ipT(rn,s = 0) = ipo(m) implies 
qris = 0) = qris = 0) = 0. The point {q, q) = (0, 0) is a stationary point of H for all values of s, hence for all (finite 
when TV — >■ oo) values of the annealing time T the solution of (|104D is qxis) = 9t(s) = 0. In particular when s = 1 
the final energy is efin{T) = —qxiiy = 0. 

2. Semi- classical dynamics for the annealing towards the paramagnet 

The semi-classical analysis of the annealing towards the paramagnet is more conveniently performed in the x- 
diagonal basis. We write \(J)t{u)) — ^ (/)t('7z, u)|to; 0)2, with (l)T{m,u) = e^^'PT(rn,u) ^ ^^^ obtain from the 

meM" 



Schrodinger equation ((7^ that (^t evolves according to 



i dipT{m,u) 



T 



du 



= -um - (1 - u)(l - m2)P/2 (cosh(2(^^(m, u))Y , 
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(107) 



the dynamical analog of Eq. (|41l) . The initial condition corresponds to the groundstate of the — (m^)^ term, and is 
thus given in this basis by (/3t(to,w = 0) = ipQ{m). The reduction of the partial differential equation (|107l) to an 
Hamiltonian system on {(7t(u),5t(w)} follows the same lines as in the annealing towards the ferromagnet, and yields 
(see Appendix [D] for the derivation) : 



^^9t(w) = ^H(<ZT(w),5T(^),^i) = 2p(l-7.)(l-9T(^)')''^%in(2$T(^))(cos(2qT(ti)))^"' 
^^gT(w) = --^n{qT{u),qT{u),u) = u - p{l - u)qT{u) (l - qT{u)^y {cos2qT{u)f . 

where the classical Hamiltonian is 



(108) 



n{q,q,u) = -uq - {1 - u){l - q^f/^ C0s{2q)P . (109) 

The initial condition is (7t(0) = $t(0) = 0, and the final energy is computed at the end of the evolution as 

It is easy to integrate numerically the two coupled ordinary differential equations (jlOSp . and we present on Fig. 1201 
some results obtained in this way. The plot on the left panel shows the instantaneous energy density 'H{qT{u)^ qxiu), u) 
as a function of the interpolation parameter u, for several (rather small) values of the annealing time T; the agreement 
with the integration of Schrodinger equation with TV = 80 is already excellent. On the right panel we concentrate on 
the final energy density, computed from the value of the solution of Hamilton equations in u = 1, as a function of T. 
The finite size effects on the results of Schrodinger equation get stronger for larger values of T, yet their extrapolation 
with a correction term of order 1/N is in very good agreement with the classical dynamics prediction. 




- eH,(T,iV) 




FIG. 20: Annealing towards the paramagnet of the p = 3 model, in the regime of constant times. Left: evolution of the 
instantaneous energy as a function of u, in the N ^ oo limit with various values of T (independent of A'"). The lines are the 
results of the integration of Hamilton equations of motion. The symbols are obtained via the integration of Schrodinger equation 



with N = 80, for the T values on which they fall on. The T ^ oo line is the ferromagnetic energy for u < its 



and its 



continuation with the iso-density argument for larger values of u. Right: the final energy at u = 1, as a function of T. The solid 
line has been obtained via the integration of Hamilton equations of motion, the other lines are the results of the Schrodinger 
evolution for various values of N. The A — > oo extrapolation was made with fits of the form efln(T, A) = efin(r) + x{T)/N. 
The horizontal dashed line is the asymptotic value enn for the T ^ oo limit (taken after A'^ — >■ oo), discussed in more details 
in Sec. IIVE3I The inset shows a zoom on the small T regime, for which the finite size effects are very small: the data for 
A = 100 are indistinguishable from the results of the Hamiltonian formalism. 



3. The long time limit of the annealing towards the paramagnet 



Let us now discuss the behavior of the final energy density for large values of T (yet finite with respect to N) . We 
expect that this large T limit matches the small t limit of the exponentially large time regime studied in Sec. IIV D 31 
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FIG. 21: Phase portraits of the classical Hamiltonian (|109[) for p = 3; from left to right u — 0.3, u = 0.57, u = Usp = 0.6, 
M = 0.7. 



namely that 



lim efin(T) = lim efi„(r) 



efin 



(110) 



In other words we do not foresee an intermediate scaling regime, as far as the energy density is concerned, between 
the constant times and the exponentially large times regimes. 

The intuitive explanation of this statement, in terms of the gap structure in the spectrum of the quantum Hamil- 
tonian, is the following. The gaps encountered on the metastable continuation of the ferromagnetic groundstate are 
exponentially small until the spinodal Ugp = 1 ~ Ssp is reached, thus for any finite T no turning on the crossing para- 
magnetic states can be performed before Usp. Around Ugp there are some polynomially small gaps that would need a 
polynomially growing time T to be resolved. However these gaps do not extend to values of u strictly greater than 
Msp, in the thermodynamic limit. Hence in the limit of large T, taken after the thermodynamic limit, the evolution 
should follow the paramagnetic energy levels that join the spinodal point, and hence lead to a final energy density 

efin- 

We shall give now a more quantitative justification of the statement (|110p . and characterize the asymptotic correc- 
tions efin (T) — efin as T — 7> oo, by analyzing the classical mechanics problem defined in Eqs. (|108ll09p . The large T limit 
of these equations corresponds to an adiabatic classical mechanics evolution, and we shall thus use the tools from the 
theory of classical adiabatic invariants [7l|, [72] . Consider first the phase portraits of the classical Hamiltonian (I109P , 
plotted on Fig. [^H For u < Ugp the classical Hamiltonian has a local minimum in {q,q) — {q^,{u),0), where q^,{u) is 
given in terms of the longitudinal magnetization m.^, of the ferromagnetic state by q^,{u) = iy]--m7(s = l^^ip {q* 
is in fact the associated transverse magnetization). The corresponding value of H is efm(s = 1 — m). In consequence 
the classical mechanics evolution has closed trajectories around this minimum, as can be seen on the first two panels 
of Fig. 1211 The initial condition <7t(0) = 5t(0) = corresponds to this minimum in u = 0, hence for T — >■ cx) the 
evolution follows this moving minimum (with corrections of order 1/T that will be discussed below), and reaches the 

point of coordinates ((jsp , 0) at Us 



-''sp 



(we denote qsp = 



1 - m% = 1/Vp^T) 



At the spinodal reached in Usp the 



ferromagnetic metastable state disappears; in this context this translates into the absence of such closed trajectories 
for u > Usp (note however that the Hamiltonian is 7r-periodic in q) , see the two last panels of Fig. I21I The T — > cxd 
evolution for u > Usp can be understood in terms of classical adiabatic invariants. Let us recall that these are quan- 
tities that depend on (g, q, u) and that have small variations along a trajectory solution of Hamilton equations, in 
the limit where the Hamiltonian of the system has a slow explicit time-dependence with respect to the instantaneous 
motion of the system, i.e. here in the large T limit. The simplest adiabatic invariant (conserved with corrections of 
order T~^) corresponds to a passage to action-angle variables, and reads 



I{q,q,u) 



'dq' , 



(111) 



■H{q' ,q' ,li)=H{q,qM) 



where the integral is performed over a trajectory starting in (g, q) that corresponds to Hamiltonian conservative 
evolution for a fixed value of u. Note that in this case the adiabatic invariant depends on (q, g, u) only through (e, u) 
where e = H^q, q, u) is the fixed energy on the trajectory. For the lines of the phase portraits that reach the points 
q = ±7r/2, this quantity can be computed by expressing g' as a function of q' and e — 'H{q' , q'). Inverting the relation 
(|109p one obtains 



I{e,u) = 2 




uq 



(l-u)(l-g2)p/2 



dg' 



(112) 
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FIG. 22: Left; the solutions qriu) of Ea. (|108|) for p = 3, T = 50, 100,400, with their adiabatic limit q,{u) when u<Uap. Right: 
similar datas for T = 100, 400, 800 plotted with the rescaling defined in Eq. pi9[) . together with the tritronquee solution of the 
Painleve equation. 

where gmin and gmax denote the extremal points of the trajectory. A moment of thought reveals that for u > Ugp this 
quantity is proportional to the integrated density of states 'Dq: compare it with the expression of Vq in (pi)) . and 
the semi-classical solution of the eigenvalue equation expressed in the x-basis given in Eq. (j42[) . This shows that the 
conservation of adiabatic invariants in the T — > oo limit is strictly equivalent to the iso-density argument used in the 
analysis of exponentially large timescales. Hence the classical mechanics adiabatic evolution between u — Ugp and 
u = 1 brings the system to the final energy egn, defined by I?o(ssp, esp) = 2?o(0, efin)- 

We shall now discuss the behavior of efin(T) — efin in the large T limit. One can expect some generic corrections 
of order T^^ to arise because of the imperfect conservation of the adiabatic invariant for large but finite T. However 
these effects are subdominant with respect to the singular corrections due to the bifurcation transition of the classical 
mechanics system at Wgp uM- The left panel of Fig. [52] displays the functions qriu) for several values of T. For large 
enough T they indeed follow with a good approximation q■^,{u) for u < Ugp, while they have an oscillating behavior for 
u > Wgp, in agreement with the shape of the phase portraits. The critical regime around Ugp plays however a crucial 
role, as will be explained now by reconsidering in a quantitative way the reasoning above. 

For u < Usp, the expansion of the Hamiltonian (jl09p around its minimum in {q,q) = {q^{u),0) yields 



'H{q,q,u) = efi„(s = 1 



u) + ^g{u)uj{uf{q- 



q*{u))^ 



2.9 («: 



~2 

■q 



0{{q 



q*{u)f,t, 



{q-q-MW) 



which corresponds to an harmonic oscillator centered in q^ (u) , with mass and pulsation given by 



5(w) = 



4p(l-M)(l-(7,(w)2)P/2 



uj{u) = 2p{l - u){l - q4u)^)2-^^l-{p-l)q,{uy 



(113) 



(114) 



In the large T limit one can set up an expansion for the trajectory of an harmonic oscillator with slowly varing 
parameters (here q^,,g,ui), under the form of oscillating terms of pulsation Tuj{u) (in the slow time u) multiplied by 
slowly varying terms. At the leading order, and taking into account the initial condition (7t(0) — 9t(0) = 0, one 
obtains 



qAu) ^ q.^{u)-^q'MJ jf\ 

1 y uj(0)g{u)uj{u) 



sin T / du'uj{u')] +0{T~^) , 



Mu) ^ \9{uHSu) - ^^:(0) J ^^^^^j;^^"^ COS (tJ\u'. in') ] + Q(T--) . 



(115) 



(116) 



This expansion is only valid for u < Ugp, because uj{u) vanishes as u 



In this limit the harmonic potential is 



not confining anymore. To continue the description of the evolution towards larger values of u we shall now expand 
the Hamiltonian around its bifurcation, and look for a scaling function that will describe the neighborhood of the 
singularity. We write: 



'H{q,q,u) = esp-efnj(ssp)(u-Usp) + ^ q^ - -^ttpiq - qspf - bp{q - qsp){u ~ Usp) + 



(117) 
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where ap and bp are two positive constants depending on p that can be obtained as partial derivatives of H in 
(qsp, 0, Msp)- The mechanical interpretation of the last three terms is a particle of mass g{usp), evolving in an energy 
potential that has a constant cubic term and a linear term whose sign changes as u crosses Usp, thus provoking the 
disappearance of its stable minimum. Eliminating q from the Hamilton equations of motion that follows from this 
truncated expansion leads to 

^^(Mu) - g.p) = ^^(9t(-) - 'Zsp)^ + ^(« - ^sp) . (118) 

We define a scaling function y{t) with 






3/5 

1 

1/5 (l)*P*^^^-^ (119) 

_ (p-i) 10 ^ ' 






The scalings with T of these changes of variables are chosen in such a way that the three terms of Eq. (|118l) are of 
the same order. The constants Cp and Cp are more arbitrary, and have been chosen here in order for the scaling 
function y{t) to be solution of the canonical form of the first Painleve equation, y"{t) = 6y{t)^ + t. We have only 
given Cp explicitly above as Cp will not appear in the final result. There exists of course an infinite family of 
solutions of the Painleve equation, selected for instance by the value of {y,y') at a given t. In our case the solution 
will be selected by a matching argument between the u — >■ Wgp limit of the first regime u < Ugp described by Eq. (IllSp , 
and the t — )■ — cx) limit of the regime described by the Painleve equation (in a neighborhood of Ugp of order T^'^'^). 
To expand Eq. (J115I) we note that in the limit u — >■ Usp one has q>^{u) = (jsp — y/2bp{usp — u)/ap + 0{{usp — u)^'^) and 
uj{u) ~ est (usp — u)^'^, where here and in the following we denote est positive constants whose precise values are not 
necessary for the reasoning. These expansions yields, in terms of the rescaled variables y and t, 

y{t) ~ -A/-^-cstr-i/2(-i)-i/8sin(r(cst-cst(-t)5/4)) . (120) 

t— >— oo V 6 



The leading term — y^—t/6 is common to several solutions of the first Painleve equation; however we note here that 
the amplitude of the oscillating term vanishes as T -> oo (for a fixed large t), hence the scaling function should be 
given by a monotonous solution with the — a/— t/6 asymptotic behaviour. This was shown in [74| to imply that y{t) 
is Boutroux [75J tritronquee solution. This solution was studied in great details in [74| . in particular the location of 
its smallest real pole to was determined numerically with great accuracy, and found to be to = 2.3841687 .... On the 
right panel of Fig. [21] we compare the tritronquee solution of the Painleve equation (determined numerically with its 
values (j/(0),2/'(0)) given in [tI]) with the curves qriu), rescaled according to (|119p . Their agreement improves as 
T increases, as expected for a scaling function. One can compute the instantaneous energy in the regime described 
by the Painleve equation, namely for u ^ Ugp + T~^^^t/Cp with t < to, and find from (|117p that it is given by 
Bsp — efjjj(ssp)('u — Usp) + 0{T~^/^). Let us also compute the integrated density of states associated to such energies, 

Vo [s = Sep - r-4/5-L,e = e,p - e',,^{s,p)T-'"'-LA ^ Vo{s,p,e,p) + MpT-^/'^ , (121) 

\ Cp Op J Op 

where Mp was defined explicitly in Eq. (|M|) . As explained above the conservation of mechanical classical invariants 
corresponds to the conservation of the integrated density of states for u > Usp- Our prediction for the large T 
behaviour of the final energy density thus reads 

efin(T) ~ efin + 2Afp^r-4/5 . (122) 

Cp 

Indeed the largest violation of the conservation of the adiabatic invariant is obtained by taking t -^ to, the limit of 
existence of the scaling regime described by the Painleve equation. 

It is rather peculiar that a scaling function matching two different regimes is defined only on a part of the real 
axis (here t < to)- In fact at the end of the Painleve regime the values of u are still close to the singularity 
{u — Usp = 0{T~"^^^)), hence the periods of the orbits encountered at those times are divergent. It has been shown 



37 



I 




ln(T) 

FIG. 23: The large T limit of the final energy density for the annealing towards the paramagnet of the p = 3 model. The 
symbols have been obtained via the integration of Hamilton equation of motion, for T as large as 24000. The two lines are of 
the form efln(T) — ea-n = aT~^'^ + bT~^'^ + cT~^. In both cases a was fixed by the analytical prediction from Eq. (|122|) . the 
function /i was obtained with c = and using & as a fitting parameter, while for /2 we fitted the data with both b and c. 

in [73| how to deal with this third regime of time, that matches the t —^ to limit with u = Wgp + £, where e is arbitrary 
small but independent on T. In particular it was found that the additional corrections to the adiabatic invariant due 
to this regime are of order T~^/^, i.e. asymptotically neglectible with respect to those we have computed, yet larger 
than the regular T~^ corrections to the action adiabatic invariant. 

We have checked the analytical prediction (I122p against numerical integrations of the Hamilton equations of motion, 
and present these results on Fig. [23l We could not achieve a good agreement with the data using only the form (112211 : 
indeed, even for the largest times T = 24000 we could reach, the subdominant correction term of order T~^/^ [73| 
is comparable to the leading one (the difference between the two exponents 4/5 and 5/6 is tiny). Including this 
correction term as a fitting parameter yields a very good agreement with the data, that is further improved with the 
inclusion of the regular T~^ corrections. We have also checked for other values of p a similar agreement with the 
prediction of Eq. ([T^ . 



Even values of ; 



Let us finally discuss the annealing for even p models that was left aside in the previous discussion. As explained 
at the end of Sec. IIII Al the models with even p enjoy an additional symmetry, the conservation of the parity of the 
magnetization in the x basis. This implies that the dynamics of the even p> A models has exactly the same properties 
as the odd p > 3 cases. Indeed the dynamics is confined to the subspace of parity equal to the one of the groundstate. 
In that subspace the ferromagnetic levels are unique, and all the structure of the gaps in the spectrum is qualitatively 
the same as for odd p > 3 models. 

The case p = 2, studied in j37l - l4C)l |. is on the contrary very different. Consider first the annealing towards the 
paramagnet. The only relevant timescale, as far as the energy density is concerned, is the one of finite T when 
TV — )■ 00. Indeed, as explained at the beginning of Sec. IIV E 3[ resolving the gaps of order N~^^^ encountered 
around Uc = 2/3 (hence considering interpolation times of order iV^") is necessary only to end up the evolution in 
the groundstate, not to reach energy densities equal to the one of the groundstate. On the finite T timescale the 
semi-classical analysis of Sec. IIV E 21 is thus relevant, and allows to cover the full range of energy densities between 
efin{T = 0) = and efjn{T -^ oo) = — 1. Moreover the large T corrections to the energy density are much less singular 
than for p > 3, because the bifurcation at u^ is of a different type. The scaling regime is described by the second 
Painleve equation [3^,|40] instead of the first one, and this should lead to corrections of the form enniT) ~ — l + cst/T. 

The annealing towards the p = 2 ferromagnet has a much richer structure. For all finite T, in the thermodynamic 
limit, the semi-classical analysis of Sec lTVE ll predicts that efin{T) = 0. Indeed the initial condition {qris = 0),q{s = 
0)) = (0,0) is a stationary point of H for all values of s, even though it becomes unstable for s > Sc- Reaching 
non-trivial final energy densities thus requires interpolation time-scales that grows with N; their precise scaling is a 
delicate problem that we leave for future work. Indeed a preliminary treatment, within the formalism of this paper, 
reveals that the gaps that close along the line e — — (1 — s) (corresponding thermodynamically to the unstable solution 
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m = of Eq. (|lip) do so as 1/lniV, and not polynomially in N as happens for the groundstate. Such a logarithmic 
behaviour was aheady discussed in [4J, [7a] ■ 

V. CONCLUSIONS 

Let us give a partial summary of this work and propose a few directions for future research. One of our main results 
in the statics part of the paper is the formula (|62p that gives the exponential rate of closing of the gap at a first-order 
phase transition (i.e. for p > 3 in this class of models), under the form of a semi-classical tunneling amplitude between 
the paramagnetic and ferromagnetic states that cross at the transition. At a second-order phase transition (here for 
p = 2) we recover the results of 30;, i57|] on the polynomial scaling N~^/^ of the gap, from a matching between the 
square-root closing of the finite gap in the paramgnetic phase (see Sec. IIIIF|) and the behavior of the exponential 
splitting of the two ferromagnetic groundstates around the transition (studied in Sec. lIII(T2|) . 

The detailed description of the spectral properties of the models studied here, in particular the density of states and 
the rate of closing of exponentially small gaps, relies on the analysis of the solutions of the semi-classical eigenvalue 
equation (j28p . It is actually straightforward to write its generalization, and thus to perform the same subsequent 
steps of analysis, for any model of spins whose Hamiltonian only depend on the total magnetizations m^, m^, m^. For 
concreteness let us give these generalizations for three examples: 

• In the LMG model the Hamiltonian reads H/N — —Tfh^ — ^x{rh^Y ~ lyi^^Y ^ ^-^d the generalization of (pS)) is 

e = -rTO-7^(l-m2)cosh^(2(^'(TO)) + 7;^(l-m2)sinh^(2(y9'(m)) . (123) 

For the density of states this should yield formulas equivalent to those obtained in [3^ . 

• For the models of [35] , where the interactions along two axis are raised to arbitrary powers, one can write the 
Hamiltonian as H/N = — 7z(m^)'' — 7x(to^)'' and the eigenvalue equation, in the thermodynamic limit, as 

e = -73 mP - 7:r(l - to^)^ cosh(2^'(TO))P' . (124) 

• The authors of [36| introduced an antiferromagnetic coupling in the interpolating Hamiltonian of the annealing, 
under the form BjN = -s[A(m^)P - (1 - A)(m^)2] - (1 - s)m^ . This yields 



e = -sXmP + s(l - A)(l - m^) cosh{2ifi' {m)f - (1 - s) y/ 1 - m'^ cosh{2ip' (m)) . (125) 

In the Section ITVl devoted to the annealing dynamics of the fully- connected p-spin models we have analyzed the final 
energy density Cfin after an evolution on a time T. Our analytical results have been obtained in the thermodynamic 
limit; let us emphasize the necessity, in this limit, to define precisely the scaling of T with the system size N. The 
results, and the methods employed to derive them, are indeed very different according to the timescale investigated. 
In Sec. IIVDI we studied annealing times T growing exponentially with N, in terms of the Landau-Zener mechanism 
controlled by the exponentially small gaps encountered by metastable states. The regime where T is kept fixed while 
the limit A'^ — > oo is performed first was analyzed in Sec. lIVEi via a reduction to a classical mechanics problem [43|. 
We have argued that these two regimes are the only relevant ones for p > 3, and as far as the energy density is 
concerned; resolving finite (extensive) energy differences would require in some cases the study of an intermediate 
timescale, with T growing polynomially or logarithmically with N. An outcome of our analysis is the crucial role 
played by spinodals in the annealing of mean-field models encountering a first-order transition: in the limit where T 
is large but finite with respect to N, or exponential with A^ but with an infinitesimal growth rate r, an annealing 
follows the metastable groundstate until its disappearance at the spinodal, and reaches at the end of the evolution an 
excited energy density Cfin corresponding to the state that crosses the metastable state at the spinodal. This energy 
separates what can be achieved on sub-exponential times (e > efin), from the range of energies egg < efin < Sfm that 
require an exponentially large annealing time to be reached. In the models studied here the paramagnetic state is 
always metastable and has no spinodal, hence the annealing from the paramagnet has a trivial finite-time regime 
(sfin — 0); this motivated the complementary study of the annealing in the reverse direction (from the ferromagnet 
to the paramagnet) which exhibits a non-trivial boundary efin between the two timescales. 

As we already emphasized the models studied in this paper are only toy-models as far as the difficulty of finding 
their groundstates is concerned; from this point of view both directions of the annealing (from the paramagnet to the 
ferromagnet or viceversa) are equally relevant. We conjecture that some of the results we obtained may remain true 
for the quantum annealing of more difficult combinatorial optimization problems as those of |17l . Il8l | . In particular the 
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scalings of the final energy density for tlie small r limit of exponentially large time-scales (see Eq. I\99\i ) and for the 
large T limit of constant timescales (see Eq. (jl22p ') with the exponent 4/5 could be generic for all mean- field models 
encountering a first-order transition followed by a separate spinodal along their quantum annealing (the location of 
the spinodals for the XORSAT problem were determined in [131): be there fully-connected or diluted, as long as they 
are mean-field. Indeed in these combinatorial optimization problems the paramagnetic state from which one starts 
the annealing procedure has a spinodal limit of existence (exactly as the ferromagnetic state of the toy models studied 
in the present paper). One of the several open questions in this context would be the determination of Cfin, in other 
words the generalization of the iso-integrated density argument that applies only to the fully-connected models whose 
Hilbert space can be decomposed in disconnected spin sectors. One possible road for this calculation in the context 
of diluted mean-field models would be the quantum extension of the "state following method" [73| (related to the 
Franz-Parisi potential J78|), that answers a similar question for classical annealing dynamics. 

In the design of a quantum annealing algorithm there is some freedom in the choice of the initial Hamiltonian 
Hi (it should however have a groundstate that is easy to prepare, and its construction should not assume a detailed 
knowledge of the sought-for groundstate of the final Hamiltonian H{). To avoid the phase transitions that appear when 
Hi is a transverse field it was for instance proposed in 79] to randomize the direction of the transverse fields on each 
spin. Very recently another proposal was to include antiferromagnetic couplings in the interpolating Hamiltonian 13611 : 
in this way it is possible to avoid the first-order phase transition by making a detour in the (s, A) plane (see also [32| 
for a similar phenomenon). The annealing towards the ferromagnet for p > 3 studied in the paper was particularly 
inefficient because the groundstate of the initial Hamiltonian (the transverse field —ffv^) remained metastable all the 
way to s = 1. One can thus wonder whether taking a ferromagnetic coupling —{ffi^Y with p' >2 (this corresponds 
to the models of [3a|) would help. The answer is no, the metastability until s = 1 persists for all values of p', as 
long as p > 3. Instead the antiferromagnetic coupling (m^)^ introduced in [3a] helps the annealing because their 
groundstate |0; O)^; has a much larger overlap with the groundstate |1; Q)z of the target Hamiltonian —{fh^Y than has 
the groundstate |1; O)^ of the transverse field. 

In the fully-connected ferromagnetic models studied in this paper the condition for a thermodynamic first-order 
transition (i.e. p > 3) coincided with the existence of exponentially small gaps at the transition. There can however 
be exceptions to this rule: the case p — p' — 2 oi [35] exhibits a discontinuity in the derivative of the groundstate 
energy but no exponentially small gaps. This peculiarity is due to the coincidence of the spinodals with the first-order 
transition, the states that cross become unstable right after the transition. A similar situation was shown to happen 
in antiferromagnetic chains of odd lenghts with periodic boundary conditions [80| . 
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Appendix A: Large p expansion of the closing rate of the gap 

This appendix is devoted to the derivation of the asymptotic expansion (1631) for the rate of the exponential closing 
of the gap at the first-order transition, in the large p limit. Let us first compute the limit of a^. Simplifying the 
expression (15^ with m^ = \^ Sc = 1/2, Cc = —1/2, one obtains: 



1 /"^ . f 1 A 1 /"^ w N ln2 
- / dmach — , = — / dm argtanhrn = 

2 70 \Vl^^) 2 7o ^ ^ ^ 2 



as argued for in [3J]. For the computation of ap at order 1/p the corrections to rric and Sc given in Eq. (jlSp are 
actually irrelevant and one has: 



"'= U' "■"-'■ (^''-""O^K?) 



2^*! Jo \s/i~^l V* A^l^^/ \V' 



(A2) 
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In general one has -4-acha:: — X) -tti '^»,fc (\_-i\i-i/2 1 ^nd therefore the fc— th integral above is found to be: 

where we used that '^^Ci^k = limaj-s-oo a;'^3prach(x) = limx^oo x'^ -j^ ^^ x = (— l)'^^^(fc — l)!x~'^. We obtain finally 
the expansion of Eq. (|53)) : 

_ ln2 1 ^ (fc-1)! /^ 

"^"2 2p^ /fc.fc! ^ U' 

''-' (A4) 

= i^_lvl /^\ ln2_^ fV 

2 2pf-^_^k^^ \p^ J 2 12p ^ \p^ 

Appendix B: Technical details on the simplified model 
1. Statics 

We justify in this appendix the formula (|85p for the rate of closing of the gaps of the simplified model, along the 
metastable continuation of its groundstate for s > 1/2. We look for an eigenvector of H,{s), with an eigenvalue Ne, 
under the form ^ 0('7i, s, e)|?n; O)^;. The coefficients of this decomposition are solutions of 

e(j){m,s,e) = —sm(t){m,s,e) — {1 — s)D,n /^ D,n'(j){m' , s,e) . 

For s = the lowest eigenstate is given exactly by (j){Tn, 0, —1) = I?™, which can be written at the leading order in the 
thermodynamic limit ^-'^voi™.)^ with (/pq given in Eq. ([Ml) . For < s < 1/2 we construct an approximation <l){m,s) 
of the groundstate eigenvector as 

^(m, s) ~ (?^(m,0,-l) 
0(m,s) = , with (/)(77i, s) = ^-— . (Bl) 

Il0(s)ll ^ ~ T^"^ 

Indeed one has : 

^^^^/.„„^_ (/.(m,0,-l)_ n V- n 0(m',O,-l) 



1 - -T-^m I 0(to, s) = ~ — = -An 2^ A 






m'GA^^ 



(B2) 



Thus 



(B3) 



- (1 - s)(/)(to, .s) = -STO0(m, s) - (1 - s)Ari ^ i:'„'(/)(m',s) + ( — 7=p ) 

This shows that, for s < 1/2, the groundstate has energy close to —(1 — s) and takes the form 

4>{m, s, -(1 - s)) = e-^^"im)+oiVN) ^g^^ 

For s > 1/2, there appears a divergence in the definition of (j){m,s) at to = ^— ^, a sign of the avoided crossing 
with an eigenvector localized near to, = —^^ in the x-basis. This divergence is lifted by constructing the symmetric 

and antisymmetric combinations of these two quasi-eigenvectors. Let (/)±(to, s) — -y= ({1 — (5,„ i^^)<j){m, s) ± S^^ i^s j . 

Then one can see that 4'± still satisfies (jB2p . and thus correspond to two quasi-eigenvectors at the location of the 
avoided crossing. At the leading exponential order the gap between the two eigenstates that crosses for some value 
of s is given by the overlap between the metastable state of eigenvector close to e~^^°^ '' and the localized state in 
TOo(s) = -!— ^. As a consequence 7»(s) = (po((l ~ s)/s), which explains the origin of Eq. 
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2. Annealing with a sub-exponential interpolation time 

In this appendix we present an analysis of the anneahng of the simplified model with an interpolation time growing 
sub-exponentially with N. The Schrodinger equation on the vector \(J)t{s)) = ^ (/)t(77i, s)|to; 0)^; reads : 

^ dM^.s) ^ _sNmMm. s) - N{1 - s)D„, Y. D^n'<l>T{m', s) = -sNmMm, s) ~ N{1 - s)/t(s) , (B5) 
where we introduced fris) = J2 Dm(f>T{m, s). By summing (jBSp over m we obtain: 

^^^^ = -sN{m)s.T-N{l-s)fT{s) (B6) 

1 as 

with {m)s.T — X] fnDm(f)T{Tn,s). Note that {ra)s=o,T = 0. Assume first that one can neglect {m)s^T in (|B6I) . 
Then using the initial value /t(0) = 1, one obtains /t(s) = e*^-^'""^ /'^\ Replacing into (jB5|) gives: 

^d^rKs) ^ _^^^^^(^^^) _ ^(^ _ ,)e.^T(.-.V2) . (B7) 

T as 

A solution of the associated homogenous equation is (pj^ {m, s) — e^"^^ Nm/2^ Writing the solution of the complete 
equation (jB7|) as (/)t(?ti, s) ~ Xrim, s)(j)rp (m, s) leads to; 

l^MUhfl = _iv(l - ,)e-^#^e^^^(-^^/^) , (B8) 

T ds 

with Xrim, s — 0) — Dm- Let us now write Xrim, s) — hj^.Ti'm, s)e^^^^™'^\ and assume that lim -i- In ST{rn,s) _ q 

A^ — ^oo 

and lim jj lnft,jv xi^i^, s) — 0. Then it is easily found that gri'm, s) — —(po{'m) — iTs^m/2 + iT{s — s^/2), and thus 

(f>T{'m,s) — (pTiTTL^s — Q)hM,T{in,s)fT[s). The two conditions above then reduces to limjv_>.oo "^InT — 0, that is, 
that one considers sub-exponential times. Finally, it is easy to check that in this regime {m)s^T — 0, and thus that 
our derivation is indeed self-consistent. 

To summarize, in this case, 0t(?7i, s) is up to subdominant corrections equal to (j)T{in, 0) times a phase independent 
on m, and the final energy efin(r) is thus identically zero. Therefore we showed that for the simplified model: 

sup lim Cfinfr ^ N'',N)^ lim efi„(r) = lim Cfinfr) = , (B9) 

a T,N—yoo T->-oo t->-0 

where the last term comes from the analysis of (small) exponential times of Sec. IIVD II 

Appendix C: Numerical integration of the Schrodinger equation 

In this appendix we explain the details of the procedure we used for the numerical treatment of the finite N dynamics, 
inspired by 81, 82]. Solving the time-dependent Schrodinger equation (I75p amounts to compute the evolution operator 
U{0, 1) = T ( e^^-^Jo i{{s)ds j .^^j^gj-g J- denotes the time-ordering operation. It is convenient numerically to break the 

time interval s e [0, 1] in n intervals of length As — 1/n, with equidistant discrete times Si = {i — l)/n. This allows 
to write: 

n n 

We are interested in the particular case of a linear dependency of H{s) on s: H{s) = {l~s)Hi+sH{. The approximation 

^ / _,T-^4^gA / -.T ^'^->t-^°^ gA (C2) 

= Uas{s) 
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gives rise to an error in operator norm |1j4|1 = sup^f ||js:=i|| ll^-'^ll bounded by |8]| : 

\\Uis,s + As)-UAsis)\\ <\\[Hi,Hf]\\^^+OiAs^) = OiNTAs^) . (C3) 

Indeed in all the cases of interest here the conrniutator of the initial and final Haniiltonian has a norm of order N. 
We define the approximate evolution operator U{0, Si) = ri7=o^'^'*('*i)- '^^^ triangle inequality 



(C4) 



\\U{0, S,+ i) - U{0, S,+ i)\\ = \\U{0, s,){U{s„ S,+ i) - Uas{s,)) + {U{{), S,) - U{0, S^))UAs{s^)\\ 

<\\U{s„S, + As)-UAs{Sr)\\+\\U{0,S,)-U{0,.S,)\\ 

leads by recurrence to 

\\U{0, l)-U{0, 1)11 < 0{nNTAs^) = 0{NT/n) . (C5) 

One can thus replace the exact evolution operator U{0, 1) by its approximation U(0, 1) with a precision of order e in 
the evaluation of intensive observables if the number of discretization steps n is of order NT/e. 

Let us evaluate the total complexity of the procedure. The dynamical evolution occurs in the fully symmetric sector 
of the Hilbert space, hence all operators are actually matrices of size N + 1. For the evolution towards the ferromagnet 
Hi = —Nfff, Hi — ~N{fh^)P, and we work in the basis where m^ is diagonal. We do not compute all the matrix 
elements of ZY(0, 1), but rather its product with the initial state |0t(O)), a column vector of size A^ + 1. For each 
time increment Si — > s^+i we have to multiply the (approximation) of \(j)T{si)) by the two matrices in (jC2p . The first 
multiplication in (|C2p is computed in a time proportional to N as H{ is a diagonal matrix. The multiplication with 
the second term is performed with 0{N'^) operations, provided m^ (whose expression in this basis is given in Eg. llSp 
is diagonalized as an initialization step (this costs 0{N^) operations). The total cost of the computation is thus 
0{e~^ N^T) + 0{N^). In the exponentially large times regime the second term becomes irrelevant; the limitations of 
this numerical method arise from the large times investigated rather than from the sizes of the matrices themselves. 
The evolution towards the paramagnet is treated similarly, the role of Hi and H{ being simply exchanged with respect 
to the previous case. The integration of the dynamics of the simplified model defined in Sec. llVCJ is slightly easier. 
One can indeed exploit the fact that ifj = J is a matrix of rank one with its non-zero eigenvalue equal to —N, thus 

e^' = l- ^ ^ > J , (C6) 

with 1 the identity matrix. This avoids the diagonalization of the matrix Hi, and in this case the total complexity 
of the integration is 0{e~^N'^T), the multiplication of a rank one matrix with a vector being computable with 0{N) 
operations. 

Appendix D: Derivation of Hamilton's equations of motion 

In this Appendix we give some details of the derivation of (|104p from (|103p . A similar computation can be found 

in m. 

The equation (jl03p on ipT{fTi, s) is of the form 

where 'H{q,q,s) is a smooth real function of its parameters. We write ipT{'m,s) — g{m,s) — i9{m,s), with g and 9 
real- valued functions (from now on we keep implicit the dependence on T). We assume that there exists a continuous 
function q{s) such that 3{q(s),s) _ g ^^^ — g(q(s),s) _^ q ^^^ ^jj times s and we let q{s) — -if-9{q{s), s). Then we get: 



dm drn^ ' ^ V / Qrn 

(D2) 



ld^9{qis),s) c^ d / .dg 89 \ dniqis),qis),s) d^giqis),s) 

— ^-; — H[m,i— \- - — ,s' — 



T drads dm \ dm dm j m=n(s) 

The derivation of the condition ^ 3^ = with respect to s yields: 

ld<z(s)_ 1 '^%mt^ _dniqis),qis),s) 



drn^ 



T ds T d'^9{q{s),s) 



(D3) 



43 



In a similar way we obtain: 



This leads to 



T dmds dm \ ' dm dm / _ , ■. 

\ / m-q(s) /p^N 

dn{q{s),q{s),s) dn{q{s),j{s),s)d^d{q{s),s) 
dq dq dm? 

1 Aq{s) ^ 1 dH{q{s),s)dq{s) ^ 1 d^e{q{s),s) ^ dH{q{s),q{s), s) 
T ds T d'^m ds T dmds dq 



(D5) 



Therefore q(s) and q{s) obey indeed Hamilton's equations of motion for the Hamiltonian T-L^qjq^s). Note that the 
reality condition on % corresponds to the Hermitianity of the quantum operator iJ, and that the derivation shows 
also how Eq. (|108l) is a consequence of Eq. (J107I) . 
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